Submitted to the Astrophysical Journal 

Preprint typeset using IATpX style emulateapj v. 08/22/09 


THE EVOLUTION OF THE GALAXY REST-FRAME ULTRAVIOLET LUMINOSITY FUNCTION OVER THE FIRST 

TWO BILLION YEARS 

Steven L. Finkelstein 1,a , Russell E. Ryan Jr. 2 , Casey Papovich 3 , Mark Dickinson 4 , Mimi Song 1 , Rachel Somerville 5 , 
Henry C. Ferguson 2 , Brett Salmon 3 , Mauro Giavalisco 6 , Anton M. Koekemoer 2 , Matthew L. N. Ashby 7 , Peter 
Behroozi 2 , Marco Castellano 8 , James S. Dunlop 9 , Sandy M. Faber 10 , Giovanni G. Fazio 7 , Adriano Fontana 8 , 
Norman A. Grogin 2 , Nimish Hathi 11 , Jason Jaacks 1 , Dale D. Kocevski 12 , Rachael Livermore 1 , Ross J. McLure 12 , 
Emiliano Merlin 8 , Bahram Mobasher 13 , Jeffrey A. Newman 14 , Marc Rafelski 15 , Vithal Tilvi 3 , S. P. Willner 7 

Submitted to the Astrophysical Journal 

ABSTRACT 

We present a robust measurement and analysis of the rest-frame ultraviolet (UV) luminosity function at z = 

4 to 8. We use deep Hubble Space Telescope imaging over the CANDELS/GOODS fields, the Hubble Ultra 
Deep Field and the Hubble Frontier Field deep parallel observations near the Abell 2744 and MACS J0416.1- 
2403 clusters. The combination of these surveys provides an effective volume of 0.6-1. 2 xlO 6 Mpc 3 over 
this epoch, allowing us to perform a robust search for bright (M uv < -21) and faint (M uv = -18) galaxies. 

We select galaxies using a well-tested photometric redshift technique with careful screening of contaminants, 
finding a sample of 7446 galaxies at 3.5 < z < 8.5, with >1000 galaxies at z ss 6 - 8. We measure both 
a stepwise luminosity function for galaxies in our redshift samples, as well as a Schechter function, using 
a Markov Chain Monte Carlo analysis to measure robust uncertainties. At the faint end our UV luminosity 
functions agree with previous studies, yet we find a higher abundance of UV-bright galaxies at z > 6. Our best- 
fit value of the characteristic magnitude M* is consistent with -21 at z > 5, different than that inferred based on 
previous trends at lower redshift. At z = 8, a single power-law provides an equally good fit to the UV luminosity 
function, while at z = 6 and 7, an exponential cutoff at the bright-end is moderately preferred. We compare 
our luminosity functions to semi-analytical models, and find that the lack of evolution in M* is consistent with 
models where the impact of dust attenuation on the bright-end of the luminosity function decreases at higher 
redshift, though a decreasing impact of feedback may also be possible. We measure the evolution of the cosmic 
star- formation rate (SFR) density by integrating our observed luminosity functions to Myv =-17, correcting for 
dust attenuation, and find that the SFR density declines proportionally to (l+z) -4 -^ 0 - 5 atz > 4, consistent with 
observations at z > 9. Our observed luminosity functions are consistent with a reionization history that starts 
at z > 10, completes at z > 6, and reaches a midpoint ( Xhii — 0.5) at 6.7 < z < 9.4. Finally, using a constant 
cumulative number density selection and an empirically derived rising star-formation history, our observations 
predict that the abundance of bright z = 9 galaxies is likely higher than previous constraints, though consistent 
with recent estimates of bright z ~ 10 galaxies. 

Subject headings: early universe — galaxies: evolution — galaxies: formation — galaxies: high-redshift — 
ultraviolet: galaxies 
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1. INTRODUCTION 

The past half-decade has seen a remarkable increase in 
our understanding of galaxy evolution over the first billion 
years after the Big Bang, primarily due to the updated near- 
infrared capabilities of the Hubble Space Telescope. Robust 
galaxy samples at z > 6 now include more than 1000 objects 
(e.g., Bouwens et al. 2010a; Oesch et al. 2010b; Finkelstein 
et al. 2010; McLure et al. 2010; Bunker et al. 2010; Finkel- 
stein et al. 2012b; Yan et al. 2012; Oesch et al. 2012) with a 
few candidate galaxies having likely redshifts as high as 10 
(e.g., Bouwens et al. 2011a; Ellis et al. 2013; McLure et al. 
2013; Oesch et al. 2013; Coe et al. 2013; Oesch et al. 2014; 
Bouwens et al. 2014). These galaxies are selected photomet- 
rically, primarily based on a sharp break at rest- frame 1216 
A due to absorption by intervening neutral hydrogen in the 
intergalactic medium (IGM). 

Studies of galaxies at z > 6 have revealed a number of inter- 
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esting results. Galaxies at 6 < z < 8 appear to have bluer rest- 
frame ultraviolet (UV) colors than at lower redshift, likely due 
to a decrease in dust attenuation, although the brightest/most 
massive galaxies do appear to have comparable dust attenua- 
tion at z = 4-7 (e.g., Finkelstein et al. 2010; Bouwens et al. 
2010b; Finkelstein et al. 2012b; Dunlop et al. 2013; Bouwens 
et al. 2013). Lower mass galaxies have colors consistent 
with stellar populations harboring significant metal content 
(though likely sub-Solar), thus the currently detectable pop- 
ulations of galaxies are not dominated by the primordial first 
generation of stars (e.g., Finkelstein et al. 2012b; Dunlop et al. 
2012, 2013). The structures of these galaxies are resolvable, 
though they show small sizes with half-light radii < 1 kpc, 
consistent with the evolution previously detected at lower red- 
shifts (e.g., Ferguson et al. 2004; Oesch et al. 2010a; Ono et al. 
2013). Finally, the abundance of these galaxies appears to ac- 
count for the necessary photons to completely reionize the in- 
tergalactic medium (IGM) by z ~ 6, and perhaps as high as z = 
7-8 if one assumes that galaxies at least 5 magnitudes below 
the detection limit of HST exist (e.g., Finkelstein et al. 2012a; 
Robertson et al. 2013). 

One of the key measurements is the galaxy rest-frame UV 
luminosity function (hereafter referred to as the luminosity 
function), as it is one of the most useful tools to study the evo- 
lution of a galaxy population. This measure encapsulates the 
relative abundances of galaxies over a wide dynamic range 
in luminosity. As the UV light probes recent star-formation 
activity, the integral of the rest-UV luminosity function pro- 
vides an estimate of the cosmic star-formation rate density 
(e.g., Madau et al. 1996; Bouwens et al. 2012; Madau & Dick- 
inson 2014), although this measurement is reliant on dust 
corrections. The luminosity function is typically parameter- 
ized with a Schechter (1976) function with a power-law slope 
at faint luminosities, and an exponentially declining form at 
the bright end. Comparing the shape of the luminosity func- 
tion to the underlying dark-matter halo mass function, pre- 
vious studies have found that the luminosity function at z < 
6, when normalized to the halo mass function at the charac- 
teristic magnitude lies below the halo mass function at 
both bright and faint luminosities. This is generally assumed 
to be due to feedback: dominated by accreting supermassive 
black holes at the bright end (active galactic nuclei; AGN), 
and by supernova or radiative-driven winds at the faint-end 
(e.g., Somerville et al. 2008). Dust extinction can also play a 
role, particularly if the level of attenuation is dependent on a 
galaxies stellar mass or UV luminosity (e.g., Finkelstein et al. 
2012b; Bouwens et al. 2013). Although luminous AGN are 
present at z= 6 (e.g., Fan et al. 2006), they are exceedingly 
rare, and to date only a single quasar has been observed at 
z > 7 (Mortlock et al. 2011). Therefore one may expect the 
degree of the exponential decline at the bright end to become 
weaker with increasing redshift. In addition, robustly quanti- 
fying the bright end of the luminosity function can allow us to 
gain physical insight into how these distant galaxies turn their 
gas into stars, as the star-formation timescale is a significant 
fraction of the age of the universe, therefore enough time has 
not yet elapsed for feedback to bring these galaxies into equi- 
librium. A change in the star-formation timescale is therefore 
more readily apparent in the shape of the bright end of the 
luminosity function (e.g., Somerville et al. 2012). 

Thanks to the combination of observations from GALEX 
and the Hubble Space Telescope estimates of the UV lumi- 
nosity function exist now from z < 1 (Arnouts et al. 2005; 
Cucciati et al. 2012) out z > 8 (e.g, Bouwens et al. 2007; 


McLure et al. 2009; Bouwens et al. 2011b; Oesch et al. 2012, 
2013; Lorenzoni et al. 2013). Earlier works have concluded 
that M£j Y declines from around -2 1 at z = 3 to fainter than - 20 
at z = 8, with the faint-end slope a becoming steeper over this 
same redshift range (e.g., Bouwens et al. 2007; Reddy & Stei- 
del 2009; Bouwens et al. 2011b; McLure et al. 2013; Schenker 
et al. 2013). Flowever, in order to adequately quantify the am- 
plitude and form of the bright end, large volumes need to be 
probed, as bright sources are relatively rare. This has been ac- 
complished via a combination of ground and space-based sur- 
veys at z < 6, with a variety of studies showing conclusively 
that a single power law does not fit the data, and that some 
sort of cut-off is needed at the bright end (e.g., Arnouts et al. 
2005; Bouwens et al. 2007; Reddy & Steidel 2009; McLure 
et al. 2009). Although previous luminosity functions have 
been published at z > 6, the space-based studies have been 
based on small volumes (e.g., Bouwens et al. 2011b), and 
thus, while they can somewhat constrain the faint-end slope, 
they do not have the capability to constrain the bright end. 

Recent studies are starting to make progress at the bright- 
end. Finkelstein et al. (2013), while selecting galaxies for 
spectroscopic followup in the GOODS-N field, found an over- 
abundance of bright galaxies at z = 7. Ono et al. (2012) 
found a similar result, with their discovery of the Muv =-21.8 
galaxy GN-108036 at z= 7.2 in GOODS-N. Likewise, Hathi 
et al. (2012) found two bright z > 6.5 candidate galaxies in 
a ground-based near-infrared survey of GOODS-N. Thus, it 
appears that the abundance of galaxies at the bright end of 
the luminosity function may not be decreasing towards higher 
redshift as previously thought. Although these studies were 
based in a single field, further evidence comes from Bowler 
et al. (2014), who used new deep ground-based near- infrared 
imaging from the UltraVISTA survey (McCracken et al. 2012) 
to discover 34 luminous z ~ 7 galaxy candidates over 1.65 
deg 2 . They combined these galaxies with the results from 
McLure et al. (2013), which included deep and wide HST 
imaging over 300 arcmin 2 in the GOODS-S, UDS and HUDF 
fields, to analyze the rest-frame U V luminosity function at z = 
7. They concluded that they did see evidence for a drop-off in 
the luminosity function at the bright end, however, the drop- 
off was less steep than that predicted by a Schechter func- 
tion, leading them to postulate that the z = 7 luminosity func- 
tion has the shape of a double-power law, perhaps similar to 
that of the possible form of far-infrared luminosity functions 
(Sanders et al. 2003; Casey et al. 2014a). 

In this study, we measure the rest-frame UV luminosity 
function at 4 < z < 8 with solely space-based data, using the 
largest HST project ever, the Cosmic Assembly Near-infrared 
Deep Extragalactic Legacy Survey (CANDELS; Pis Faber & 
Ferguson). The large area observed by CANDELS allows us 
to probe large volumes of the distant universe for the rare, 
bright galaxies. With these data, we investigate the form of 
the bright-end of the luminosity function and the implications 
on galaxy evolution. In addition to the deep data in the HUDF, 
we use the CANDELS data in the GOODS-S and GOODS- 
N fields, which have not only deeper near-infrared imaging, 
but also imaging in more optical and near-infrared filters than 
the other three CANDELS fields (UDS, EGS and COSMOS). 
We also include in our analysis the parallel fields from the 
first year dataset of the Hubble Frontier Fields, near the Abell 
2744 and MACS J04 16. 1-2403 galaxy clusters. The combina- 
tion of these data allows us to select a large sample of nearly 
7500 galaxies, over a wide dynamic range in UV luminosity 
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atz = 4-8 (Figure 1). 

This paper is organized as follows. In § 2 we discuss the 
imaging data used and the catalog construction, and in §3 we 
present our sample selection via photometric redshifts, and 
estimates of the contamination. In § 4 we highlight our com- 
pleteness simulations, and in § 5 we discuss the construction 
of the rest-UV luminosity function atz= 4, 5, 6, 7 and 8. In § 6 
we discuss the implications of our luminosity function results, 
while in § 7 we compare our results to semi-analytical mod- 
els. In § 8 we present our measurements of the cosmic star- 
formation rate density, and in § 9 we discuss the implications 
for galaxies at higher redshifts. Our conclusions are presented 
in § 10. Throughout this paper we assume a WMAP7 cos- 
mology (Komatsu et al. 2011), with H 0 = 70.2 km s _1 Mpc -1 , 
VLm = 0.275 and Oa = 0.725. All magnitudes given are in the 
AB system (Oke & Gunn 1983). 

2 . OBSERVATIONS AND PHOTOMETRY 
2.1. Imaging Data 

Studying galaxies in the early universe requires extremely 
deep imaging, necessitating space-based data. Additionally, 
to probe a large dynamic range in luminosities, we need to 
combine deep studies over small areas with larger-area sur- 
veys with shallower limiting magnitudes. Our study used 
imaging data from a number of surveys covering both the 
Northern and Southern fields from the Great Observatories 
Origins Deep Survey (Giavalisco et al. 2004), with both the 
Hubble Space Telescope ( HST) and the Spitzer Space Tele- 
scope. 

The deepest imaging comes from three surveys of the Hub- 
ble Ultra Deep Field (HUDF): the original HUDF survey 
which obtained optical imaging with the Advanced Camera 
for Surveys (ACS; Beckwith et al. 2006); and the more recent 
HUDF09 (PI Illingworth; e.g., Bouwens et al. 2010a; Oesch 
et al. 2010b) and UDF12 surveys (PI Ellis; Ellis et al. 2013; 
Koekemoer et al. 2013), which obtained near-infrared imag- 
ing with the Wide Field Camera 3 (WFC3). The full HST 
dataset over the HUDF comprises imaging in eight bands: 
F435W, F606W, F775W and F850LP with ACS, and F105W, 
F125W, F140W and F160W with WFC3 (hereafter referred 
to as 7?435, V 1506 , his, z 850, T 105 , -^25> JH\^ and // 160 , respec- 
tively), which cover an area of ~5 arcmin 2 . The HUDF09 
survey also obtained deep WFC3 imaging over two similarly- 
sized flanking fields, first observed with ACS in the UDF05 
survey (PI Stiavelli; Oesch et al. 2007), referred to as the 
HUDF09-01 and HUDF09-02 fields (Bouwens et al. 2011b). 
These fields each have imaging in the F 6 06 , his, z 850 , 7i 05 , ^ 125 , 
and//i6o bands. 

The majority of our galaxy sample comes from the Cos- 
mic Assembly Near-infrared Deep Extragalactic Legacy Sur- 
vey (CANDELS; Pis Faber and Ferguson; Grogin et al. 2011; 
Koekemoer et al. 2011). CANDELS is the largest HST project 
ever, comprising 902 orbits over five extragalactic deep fields, 
including the two GOODS fields (Giavalisco et al. 2004), 
which finished in August 2013. CANDELS is composed of 
a deep and a wide survey. The deep survey covers the cen- 
tral ^50% of each of the two GOODS fields, while the wide 
survey covers the remainder of the GOODS-N field, and the 
southern ~25% of the GOODS-S field to depths ~ 1 mag 
shallower than the deep survey (the wide survey also covers 
three additional fields not used in this study; see §6.4.1 and 
Figure 15). We use ACS imaging from the original GOODS 
survey in the £ 435 , V^oe, *775 andzg 5 o bands. We use the most 


recent ACS mosaics in these fields, which is version 3 in the 
GOODS-S field, which includes all ACS imaging in that field 
prior to the ACS repair on Servicing Mission 4 in 2009, and 
version 2 in the GOODS-N field, which includes all ACS 
imaging from the GOODS survey. The CANDELS imaging 
in both the deep and wide regions of both GOODS fields in- 
cludes the Tios, J 125 and //160 bands. We add to our GOODS- 
S dataset imaging over the northern ~25% of the GOODS-S 
field from the WFC3 Science Oversight Committee’s Early 
Release Science (ERS) program (PI O’Connell; Windhorst 
et al. 2011), which also includes 7i25 and H\ 60 imaging, as 
well as the F098M (hereafter referred to as 7o9s) band. Unless 
otherwise distinguished, throughout the paper we will refer to 
7 o 98 and 7 05 together as the 7-band ( both filters probe ob- 
served 1 /xm light, but the 7o98 filter is ^half of the width of 
7 io 5 , for somewhat higher spectral resolution). 

Finally, we complete our dataset with the recently obtained 
deep HST observations near the galaxy clusters Abell 2744 
and MACS J0416.1-2403 (hereafter MACS0416) from the 
Hubble Frontier Fields (HFF) program (PI Lotz). For this 
study, we use only the parallel (unlensed) fields. Both fields 
have been observed in the R435, F 60 6, hn, 7 105 , ^125, 777 i 4 o 
and//i 6 o bands. We use these data to complement our galaxy 
samples atz= 5, 6 , 7 and 8 . 

In parallel to the primary WFC3 observations, CANDELS 
obtained extremely deep imaging in the F814W band (here- 
after / 814 ) in both of the GOODS fields. As these data were 
obtained recently, they suffers from poor charge transfer effi- 
ciency. Although algorithms have been devised to correct for 
this (Anderson & Bedin 2010), as the CANDELS fields have 
imaging in both the 2775 and zg 5 o bands, we do not include the 
CANDELS 7814 photometry in the initial photometric redshift 
fitting (though we do explore its inclusion in §3.6). However, 
we did use these very deep data during our visual inspection 
step, which was highly useful atz= 8 , where truez = 8 galax- 
ies should be completely undetected in the / 814 -band. In the 
HFF parallel fields, where the /gi 4 band is the only imaging 
covering the red end of the optical, we used these data in the 
full analysis. 

The description of the CANDELS HST imaging reduction 
is available from Koekemoer et al. (2011). These reduction 
steps were also followed for the ERS, HUDF (Koekemoer 
et al. 2013) and HFF data we use here. We use imaging 
mosaics with 0.06" pixels, and make use of their associated 
weight and rms maps. The combined imaging dataset covers 
an area of 301.2 arcmin 2 , with 5a limiting magnitudes in the 
//160 band ranging from 27.4 to 29.7 mag (measured in 0.4" 
apertures). These datasets are summarized in Table 1. 

2 . 2 . Point Spread Function Matching 

The HST imaging used here spans more than a factor of 
three in wavelength, thus the differences in point-spread func- 
tion (PSF) full-width at half-maximum (FWHM) across that 
range are significant. For example, the PSF in the GOODS-S 
Deep field has a FWHM = 0.193" in the // 160 -band, but only 
0.1 19" in the 8435 -band. A point-source will thus have more 
of its flux contained within a 0.4" aperture in the , 8435 -band 
compared to the // 160 -band. As the selection of distant galax- 
ies relies very heavily on accurate colors, and we are using 
apertures of fixed sizes (determined by the detection image, 
see §2.3) to measure photometry in all bands, this changing 
PSF needs to be addressed. 

We corrected for this by matching the PSF of the HST 
imaging to the // 160 -band image (which has the largest PSF 
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TABLE 1 

Summary of Data - Limiting Magnitudes 


Field Area 

(arcmin 2 ) 

S 435 

(mag) 

K 06 

(mag) 

Z 775 

(mag) 

2814 

(mag) 

Z850 

(mag) 

h)98/105 

(mag) 

J\25 

(mag) 

JH ho 

(mag) 

Hi 60 
(mag) 

GOODS-S Deep 

61.6 

28.2 

28.6 

27.9 

28.1 

27.8 

28.2 

28.1 

— 

27.9 

GOODS-S ERS 

41.4 

28.2 

28.5 

27.9 

27.9 

27.6 

27.6 

28.0 

— 

27.8 

GOODS-S Wide 

35.6 

28.2 

28.7 

28.1 

27.9 

27.9 

27.3 

27.6 

— 

27.4 

GOODS-N Deep 

67.6 

28.1 

28.3 

27.9 

— 

27.7 

28.1 

28.3 

— 

28.1 

GOODS-N Wide 

71.7 

28.1 

28.4 

27.8 

— 

27.6 

27.3 

27.4 

— 

27.4 

HUDF Main 

5.1 

29.5 

30.0 

29.7 

— 

29.1 

29.9 

29.6 

29.6 

29.7 

HUDF PARI 

4.7 

— 

29.0 

28.8 

— 

28.5 

28.9 

29.0 

— 

28.8 

HUDF PAR2 

4.8 

— 

29.0 

28.7 

— 

28.3 

28.9 

29.2 

— 

28.9 

MACS0416 PAR 

4.4 

28.8 

28.9 

— 

29.2 

— 

29.2 

29.0 

29.0 

29.0 

Abell 2744 PAR 

4.3 

29.0 

29.1 

— 

29.2 

— 

29.1 

28.8 

28.8 

28.9 

Zeropoints 

— 

25.68 

26.51 

25.67 

25.95 

24.87 

26.27 

26.23 

26.45 

25.95 


Note. — The magnitudes quoted are 5 a limits measured in 0.4"-diameter apertures on non-PSF 
matched images. 


FWHM) in each field. We did this using the IDL deconv tool 
Lucy-Richardson deconvolution routine, in the same way as 
Finkelstein et al. (2010, 2012b). This routine requires the PSF 
for a given band as well as a reference PSF (in this case, the 
i/i6o-band), and it generates a kernel. The PSFs were gen- 
erated by stacking stars in each field in each band, where 
the stars were selected via identifying the stellar locus in a 
half-light radius versus magnitude plane. Each star was then 
visually inspected to ensure that there were no bright near- 
neighbors, and then the stars were stacked, subsampling by 
a factor of 10 to ensure an accurate centroiding of each star 
(i.e., to avoid smearing the PSF during the stacking). Using 
these PSFs, the deconvolution routine performed an iterative 
process, and relies on the user to determine the number of 
iterations. We did this by making a guess as to the correct 
number of iterations, and then changing this number until the 
stars in the PSF-matched images in a given band had curves- 
of-growth which matched the /7]6o-band curves-of-growth to 
within 1% at a radius of 0.4". The images were then con- 
volved with the final kernel to generate PSF-matched images. 

2.3. Photometry’ 

Photometry was measured on the PSF-matched dataset with 
a modified version of the Source Extractor software (v2.8.6, 
Bertin & Arnouts 1996). Our modified version adds a buffer 
between the source and the local background cell and removes 
spurious sources associated with the distant wings of bright 
objects. Catalogs were generated independently in each of 
our ten sub-fields, using Source Extractor in two-image mode, 
where the same detection image was used to measure pho- 
tometry from all available HST filters. For most of our fields, 
we used a weighted sum of the F125W and F160W images 
as the detection image, to increase our sensitivity to faint ob- 
jects. IntheHUDF main field and the MACS04 16 and A2744 
HFF parallel fields, we supplemented this catalog with cata- 
logs using 10 additional detection images, derived by stacking 
all possible combinations of adjacent WFC3 filters. In these 
three fields, a combined catalog was made up of all unique 
sources in the catalogs, using a 0.2" matching radius. This 
allowed very blue sources that may be too faint in the / /| eo 
image to be selected in the original F125W+F160W-selected 
catalog to be included. This procedure was replicated in our 
completeness simulations (§4). To derive accurate flux uncer- 
tainties, Source Extractor relies on both an accurate rms map, 
and a realistic estimate of the effective gain. The provided rms 


map has been shown to produce accurate uncertainties, and it 
has been corrected for pixel-to-pixel correlations which occur 
as a result of the drizzling process (see Guo et al. 2013). The 
effective gains were computed for each band separately as the 
1 (2.5) x the total exposure time for the ACS (WFC3) imag- 
ing. We have previously verified that the uncertainties mea- 
sured in this manner on 1 1ST imaging are accurate (F inkelstein 
et al. 2012b). The zero-points to convert the observed fluxes 
into AB magnitudes are given in Table 1, and are appropriate 
for the dates when these data were taken. 

Following our previous work (Finkelstein et al. 2010, 
2012a, b, 2013), colors were measured in small Kron aper- 
tures with the Source Extractor Kron aperture parameter 
PHOTAUTOPARAMS set to values of 1.2 and 1.7. Finkel- 
stein et al. (2012b) found that these apertures result in more 
reliable colors for faint galaxies when compared to isopho- 
tal or small circular apertures. An aperture correction to the 
total flux was derived in the //-band and was computed as 
the ratio between this small Kron aperture flux, and the de- 
fault Source Extractor MAGAUTO flux, which is computed 
with PHOT AUTOPARAMS = 2.5, 3.5. These aperture cor- 
rections were then applied to the fluxes in all filters. To see 
if our aperture corrections accurately recovered the total flux, 
we examined our completeness simulations (discussed in §4), 
and found that after applying this aperture correction recov- 
ered fluxes were typically 5% fainter in each band than their 
input fluxes. We thus increased the flux in each band by 5% 
to derive our best estimate of the total flux (with the exception 
of the HUDF main field, where the measured correction was 
2%). 

The Source Extractor catalogs from each band were com- 
bined into a master catalog for each field. At this step, the 
observed fluxes were corrected for Galactic extinction using 
the color excess E(B-V) from Schlafly & Finkbeiner (2011) 
appropriate for a given field 1 , and using the Cardelli et al. 
(1989) Milky Way reddening curve to derive the corrections 
based on each filter’s central wavelength. We used a mask im- 
age to remove objects in regions of bad data, where the mask 
was generated using a threshold value from the weight map. 
This mask primarily trims off the noisier edges of the imag- 
ing, but it also excludes the “death star” region on the WFC3 
array when the number of dithers was low (i.e., in the CAN- 
DELS Wide regions). The areas quoted in Table 1 are those of 

1 http://ned.ipac.caltech.edu/forms/calculator.html 
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Fig. 1. — The absolute magnitude distribution of all galaxies in our redshift samples. The shaded color denotes which of the sub-fields a given galaxy was 
detected in. This figure demonstrates that while the HUDF is useful for finding the faintest galaxies, the CANDELS imaging is necessary to discover much larger 
numbers, as well as to probe a large dynamic range in luminosity. 


the good regions in these masks. Objects were also removed 
from the catalog if they had a negative aperture correction, 
which applied to a very small number of sources, primarily 
restricted to areas near very bright sources. The remaining 
objects comprised our final catalog in each field. 

3. SAMPLE SELECTION 
3.1. Photometric Redshifts 

We selected our high-redshift galaxy sample via a photo- 
metric redshift fitting technique. This has the advantage in 
that it uses all of the available photometry simultaneously, 
rather than the multi-step Lyman break galaxy (LBG) method, 
which selects galaxies using two colors, and then may impose 
a set of optical non-detection criteria (e.g., Bouwens et al. 
2007, 2011b). Another advantage of photometric redshifts 
is that one obtains a redshift probability distribution function 
(PDF), which not only allows one to have a better estimate 
of the redshift uncertainty (o- typically ~0.2— 0.3 versus 0.5 
for the LBG technique), but can also be used as a tool in the 
construction of the sample itself. That being said, initial work 
comparing the differences between galaxy samples selected 
via both LBG and photometric redshift techniques found that 
the resulting sample properties are fairly similar (e.g., McLure 
et al. 2013; Schenker et al. 2013). 

Photometric redshifts for all sources in the catalogs for each 
fields were measured using the EAZY software (Brammer 
et al. 2008). The input catalog used all available HST pho- 
tometry, with the exception of the F814W imaging in the 
CANDELS fields, which was used solely for visual inspec- 
tion (see §2.1). We used an updated set of templates pro- 
vided with EAZY based on the PEGASE stellar population 
synthesis models (Fioc & Rocca-Volmerange 1997), which 
now include an increased contribution from emission lines, 
as recent evidence points to strong rest-frame optical emis- 
sion lines being ubiquitous amongst star-forming galaxy pop- 
ulations at high-redshift (e.g., Atek et al. 2011; Finkelstein 
et al. 2011; van der Wei et al. 2011; Smit et al. 2014; Stark 


et al. 2013; Finkelstein et al. 2013). We also include an ad- 
ditional template based on the rest-frame UV spectrum of the 
young, unreddened galaxy Q2343-BX418 (Erb et al. 2010), 
as it retains characteristics expected in high-redshift galax- 
ies, such as strong Lyo emission. EAZY assumes the inter- 
galactic medium (IGM) prescription ofMadau (1995). EAZY 
does have the option to include magnitude priors when fitting 
photometric redshifts, which uses the luminosity functions as 
a prior for whether a galaxy at a given apparent magnitude 
resides at a given redshift. As we show later, there is still 
non-negligible uncertainty at the bright end of the luminosity 
function, therefore we did not include these magnitude priors 
during our photometric-redshift fitting process. 

3.2. Selection Criteria 

We selected galaxy samples in five redshift bins centered 
at z ~ 4, 5, 6, 7 and 8 with A z = 1, using criteria similar to 
our previous work (Finkelstein et al. 2012b, 2013). The cos- 
mic time elapsed between our last two bins at z « 7 and z « 
8 is ~125 Myr. This time is much longer than the dynami- 
cal time of the systems we study, and thus leaves significant 
time for evolution. However, as studies of galaxy evolution 
move towards higher redshift, this will not always be the case 
(e.g., A tz= i 3 _>i 2 = 40 Myr) thus future studies with the James 
Webb Space Telescope will need to pay careful attention to the 
choice of sample redshifts when studying galaxy evolution. 

Rather than relying solely on the best-fit redshift value, 
we used the full redshift probability distribution curves P(z) 
calculated by EAZY (where P(z) oc exp(-\ 2 ), normalized to 
unity). Our selection criteria are: 

1) A > 3.5 significance detection in both the A 25 and 
Z /160 bands. A requirement of a significant detection in 
two bands removes nearly all spurious sources, as the 
chances of a noise peak occurring in two images at the 
same position are very small (§3.8.1). This requirement 
also limits our analysis to galaxies with z < 8.5, as the 
Lyman-break shifts into the A 25 band atz= 8.1. 
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Fig. 2. — Left) The distribution of photometric redshifts in our galaxy sample. We highlight the differences in numbers between sources discovered in the 
CANDELS GOODS fields (including the ERS), and those in the combined five deep fields from the HUDF09 and HFF parallel programs. The shallower yet 
much wider CANDELS imaging dominate the numbers in every redshift bin, by a factor of 10 at z = 4—5, and ~2 at z = 7-8, though the deep fields are necessary 
for constraints on the faint-end of the luminosity function. Right) A comparison between the spectroscopic redshift and our best-fit photometric redshifts for the 
171 galaxies in our sample with spectroscopic redshifts in the literature. The red circles denote galaxies with \z spe c— z p hot\ > 1 at >3 a significance. There are 
only six such galaxies, and all have spectroscopic redshifts at z > 4. 


2) The integral under the primary redshift peak must 
comprise at least 70% of the total integral of the redshift 
probability distribution function. This enforces that no 
more than 30% of the integrated redshift PDF can be in 
a secondary redshift solution. 

3) The integral under the redshift PDF in the redshift 
corresponding to a given sample (i.e., 6.5 - 7.5 for the 
z = 7 sample) must be at least 25%, which ensures that 
the redshift PDF is not too broad. 

4) The area under the curve in the redshift range of in- 
terest must be higher than the area in any other redshift 
range (i.e., for a galaxy in the z = 7 sample, the integral 
of P(6.5 < z < 7.5) must be higher than the integral in 
any other redshift bin). This criterion ensures that a 
given galaxy cannot be included in more than one red- 
shift sample. 

5) At least 50% of the redshift PDF must be above 
z sample ~ 1 (i.e., f P(z > 6) > 0.5 for z samp ie = 7), and 
the best fit redshift must be above z samp i e - 2. 

6) The x' 2 from the fit must be less than or equal to 60. 
This criterion ensures that EAZY provides a reasonable 
fit, though in practice it does not reject many sources. 

7) Magnitude in the I 1\ W j band must be > 22. This ef- 
fectively cleans many stars from our sample, but the 
limit is bright enough such that it is still more than 
two magnitudes brighter than our brightest z > 6 galaxy 
candidate. At z = 4, we do have a few sources close to 
this limit, but only two sources are brighter than H = 
22.4. This fact, coupled with the observation that the 
very few sources at H < 22 that satisfy our z = 4 se- 
lection criteria are either obvious stars, or diffraction 
spikes, implies that this criterion should not signifi- 
cantly affect our luminosity function results. 


Of these criteria, items #1 and #2 are by far the most con- 
straining, as most galaxies which meet these criteria, with 
Zbest > 3.5 make it into our sample. Items #3 and #4 are 
responsible for putting a galaxy in a given redshift sample. 
While some of the cuts above are arbitrary, these choices will 
be corrected for as we apply these same criteria to our com- 
pleteness simulations discussed in §4. In Figure 2 we compare 
the photometric redshifts for 171 galaxies in our sample to 
available spectroscopic redshifts in the literature 2 . The agree- 
ment is excellent, with ctaz/( 1 +z) = 0-03 1 (derived by taking an 
iterative 3er-clipped standard deviation), though the number 
of confirmed redshifts at z > 6.5 is small (only five galaxies). 
The number of outliers is also small, with only six out of 1 7 1 
galaxies (3.5%) having a photometric redshift differing from 
the spectroscopic redshift by Az > 1 at >3 ct significance. All 
of these six galaxies have z spec > 4, thus no galaxies in our 
sample have a catastrophically lower spectroscopic redshift. 
In comparison, defining outliers in the same way, we find that 
the published CANDELS team photometric-redshift catalog 
has 13 outliers out of 174 total spectroscopic redshifts, for a 
somewhat higher outlier fraction of 7.5% (Dahlen et al. 2013). 
Although the fraction of galaxies with confirmed redshifts is 
relatively small, the available spectroscopy confirms that our 
selection methods yield an accurate high-redshift sample. 

3.3. Visual Inspection 

As the candidate selection process is automated, for a truly 
robust galaxy sample, we required a visual inspection each of 
our ~7500 candidate high-redshift galaxies. During the visual 
inspection, we examined the following features: 

• Is the source a real galaxy? Objects were inspected to en- 
sure that they were not an artifact, examples of which are a 

- The spectroscopic redshifts come from a compilation made by N. Hathi 
(private communication) which include data from the following studies: 
Szokoly et al. (2004); Grazian et al. (2006); Vanzella et al. (2008, 2009); 
Hathi et al. (2008); Barger et al. (2008); Rhoads et al. (2009); Wuyts et al. 
(2009); Balestra et al. (2010); Ono et al. (2012); Kurk et al. (2013); Rhoads 
et al. (2013); Finkelstein et al. (2013). 
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Fig. 3 . — Left and Center: Color-color plots. Blue circles and squares denote objects accepted as J < 26, z > 6 galaxies, with squares indicating the ones with 
half-light radii <0.17 ,/ . Cyan stars denote candidates originally selected as galaxies but reclassified as stars based on their sizes and colors. Small circles denote 
known stars from the the SpeX Prism Spectral Libraries from the 3m NASA Infrared Telescope Facility, with spectral types as indicated in the legend. Right: 
Half-light radius versus magnitude for J <26 candidates. Symbols are the same as in the other panels. No compact galaxies have colors similar to known stars 
in both color-color plots. Similar plots were used to exclude stellar contaminants at z = 4 and 5. 


part of a diffraction spike (which frequently appear in differ- 
ent places in the ACS and WFC3 imaging due to different roll 
angles during the respective observations), oversplit regions 
of bright galaxies, or noise near the edge of the images. 

• Is the aperture drawn correctly? While the small Kron 
apertures yield the most reliable colors, they are also sus- 
ceptible to “stretching” (i.e., becoming highly elongated) in 
regions of high noise or near very bright objects. For each 
source, we compared the ratio of the flux between the Kron 
aperture and a 0.4"-diameter circular aperture to that same 
quantity for objects of a similar magnitude from the full pho- 
tometry catalog. If an object had a value >30% higher than 
similarly bright sources in the full photometry catalog and the 
aperture looks to have been affected by noise/bright sources, 
we adjusted the photometry of the object in question accord- 
ingly, using the 0.4"-to-total correction of similarly-bright ob- 
jects in the catalog. In practice, these issues affected <10% 
of galaxies in our high-redshift sample. 

• Is there significant optical flux that did not get measured 
correctly? Primarily due to the issues with inaccurate aper- 
tures discussed in the above bullet, a very small number of 
sources appeared to have optical flux when visually inspected 
that was not measured to be significant in our catalog (i.e., in 
the case of a too-large aperture, the flux is concentrated in a 
small number of pixels, while the flux error comes from the 
full aperture, so the signal-to-noise is low). In these cases, ob- 
jects were removed from our sample. This step is somewhat 
qualitative, as there are cases of objects where the aperture 
appears correct, yet there is still a ~l-2cr detection in a single 
optical band. In the majority of these cases, as we are confi- 
dent in our photometric redshift analysis, we left these objects 
in the sample. During this step, we also examined 7 gi 4 pho- 
tometry for each source in the CANDELS fields, which pri- 
marily benefits the selection of z = 8 galaxies, which should 
not be visible at this wavelength. Three z = 8 candidates with 
observable /si 4 flux were removed from our sample. 

3.4. Stellar Contamination 

The most crucial step in our visual inspection is the classifi- 
cation and removal of stellar sources, as stellar contamination 
would dominate the bright-end of the luminosity function if 
these contaminants were not considered. In particular, M-type 


stars as well as L and T brown dwarf stars can have similar 
colors (including optical non-detections) as our high-redshift 
galaxies of interest, particularly at z > 6 . While some studies 
use dwarf star colors during their selection (e.g., Bowler et al. 
2012; McLure et al. 2013; Bowler et al. 2014), many use pri- 
marily the Source Extractor “stellarity” parameter to diagnose 
whether a compact object is a star or a galaxy (e.g., Bouwens 
et al. 2011b, 2014). However, the stellarity parameter loses 
its ability to discern between a point source and a resolved 
source for faint objects. To test this further, we examined the 
stellarity of sources in the CANDELS GOODS catalogs. At 
very bright magnitudes (J 125 < 24), there is a clear separa- 
tion between stars and galaxies, with objects either having a 
stellarity near unity (i.e., stars), or having stellarity near zero 
(i.e., galaxies). However, this separation becomes less clear 
at J \ 25 > 25, where the stellar and galaxy sequences begin to 
blend together. Therefore, stellarity can be an unreliable star- 
galaxy separator at J \ 25 > 25, which is similar to the bright- 
ness of our brightest z > 7 galaxies. 

While the GOODS fields cover relatively small regions on 
the sky, the potential number of brown dwarf contaminants, 
even at J \ 25 > 25, is significant. The Galactic structure model 
of Ryan & Reid (in prep) predicts the surface density of brown 
dwarfs in our covered fields. In the GOODS-S region, using 
the area covered by the CANDELS, ERS and HUDF09 ob- 
servations, we would expect ~6 stars of spectral type M6-T9 
with ./] 25 -band magnitudes between 25 and 27. The surface 
density of M6-T9 stars in GOODS-N is similar, with an ex- 
pected number of stars in the field of ~5. Thus, the expected 
number of 25 < J 125 < 27 stars of spectral type M6-T9 in our 
whole surveyed region is ~11. While this number is small, 
the numbers of brown dwarfs are expected to fall off toward 
fainter magnitudes, thus the majority of these are likely have 
J 125 close to 25. This magnitude is similar to those of the 
brightest galaxies in our sample, which dominate the shape of 
the bright end of the galaxy luminosity function. As stellarity 
is an unreliable method of identifying these sources, we must 
find an alternative method. 

Although brown dwarfs can have similar colors to z > 6 
galaxies, and can be included in the initial sample, they fall on 
well-defined color sequences, and can thus be distinguished 
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from true galaxies. Figure 3 shows two color-color plots 3 
which we used in tandem with the size information, exam- 
ining not only stellarity, but also the FWHM and half-light ra- 
dius as measured by Source Extractor, to identify stars lurking 
our sample (similar plots were used at z = 4 and 5). If a galaxy 
appeared un-resolved (defined as having a stellarity > 0 . 8 , or 
a half-light radius and/or FWHM similar to that of stars in the 
field) then we examined that object in the color-color plots 
as shown in Figure 3. If the object also had colors similar to 
a dwarf star, then we removed it from our sample. Over all 
of our fields, we had a total of 23 objects flagged as stars in 
our z > 6 samples; 18 from our initial z ~ 6 galaxy sample, 
and 5 from our initial z ~ 7 galaxy sample. These objects 
were removed from our sample. One of these stars removed 
from our z ~ 7 sample was previously flagged as a probable 
T-dwarf by Castellano et al. (2010). We examined the sub- 
set of eight of these stars which were detected in the FourStar 
Galaxy Evolution (zFourGE) medium band imaging survey of 
a portion of GOODS-S, and found that all eight have z—J 1 and 
J 1 —J3 colors consistent with brown dwarf stars (Tilvi et al. 
2013). Of these, six stars have J\ 25 > 25, meaning that our 
high-redshift galaxy selection criteria also originally selected 
~ half of the expected number of faint brown dwarfs in this 
field. Four of these six stars have Source Extractor stellarity 
measurements < 0 . 8 , thus a stellarity-only rejection method 
would have failed to remove them. We conclude that our vi- 
sual inspection step efficiently removed stellar contaminants 
from our sample, but we emphasize that the color examina- 
tion portion was crucial to exclude the faintest stars from our 
sample. 

3 . 5 . Active Galactic Nuclei 

We screened for the presence of bright active galactic nu- 
clei (AGN) in our sample by searching for counterparts in 
Chandra X-ray Observatory’ point source catalogs. In the 
GOODS-S field, we used the 4 Msec Chandra Deep Field - 
South (CDF-S) catalog of Xue et al. (2011), and in GOODS- 
N, we used the 2Msec Chandra Deep Field - North catalog 
of Alexander et al. (2003). These catalogs have average posi- 
tional accuracies of 0.42" and 0.3", respectively. To be con- 
servative, we searched for matches in each catalog out to a ra- 
dius of 1". We then visually inspected each of the 34 galaxies 
in our sample with a match. Seven objects, all with Chandra 
catalog separations > 0 . 6 ", had nearby counterparts with posi- 
tions consistent with the Chandra catalog, and thus are likely 
providing the X-ray emission; none of these sources had spec- 
troscopic redshifts in the CDF-S catalog. The remaining 27 
sources, all with separations <0.6 ", had Chandra positions 
consistent with the X-ray emission coming from the galaxies 
in our sample. Secure spectroscopic redshifts were available 
for four of these 27 galaxies in the CDF-S catalog, ofz = 3.06, 
3.66, 3.70 and 4.76. These 27 galaxies (25 fromourz= 4 sam- 
ple, and two from our z = 5 sample) were removed from our 
galaxy sample. This removal is conservative, as although the 
X-ray detections imply the presence of an AGN, it does not 
prove that the AGN dominates the UV luminosity. 

3 . 6 . Photometric Redshifts with Spitzer/IRAC Photometry 

As we will discuss below, one of the main results of this 
work is an apparent constant value of M( JV at z > 5, brighter 

3 This research has benefitted from the SpeX Prism 
Spectral Libraries, maintained by Adam Burgasser at 
http://pono.ucsd.edu/~adam/browndwarfs/spexprism 


than many previous works. It is thus imperative that we have 
high confidence that our bright galaxies are all in fact at high- 
redshift, and not lower-redshift contaminants. To provide a 
further check on our bright sources, we re-examined the pho- 
tometric redshifts of our bright galaxies with the addition of 
Spitzer Space Telescope Infrared Array Camera (IRAC; Fazio 
et al. 2004) imaging over our fields. This imaging probes the 
rest-frame optical at these wavelengths, and thus provides sig- 
nificant constraining power because the most likely contami- 
nants are red, lower-redshift galaxies, which would have very 
different fluxes in the mid-infrared than true high-redshift 
galaxies. We examined sources with M \ 500 < -21, which is 
approximately the value of Ml ) v at these redshifts, and pro- 
vides samples of 164, 85, 29, 18 and 3 bright galaxies at z = 
4, 5, 6 , 7 and 8 , respectively. 

During the cryogenic mission, the GOODS fields were ob- 
served by the GOODS team (Dickinson et al., in prep) at 3.6, 
4.5, 5.8, and 8.0 pm. Later, during Cycle 6 of the warm mis- 
sion, broader regions encompassing the GOODS footprints 
were covered by the Spitzer Extended Deep Survey (SEDS; 
Ashby et al 2013) to 3cr depths of 26 AB mag at both 3.6 and 
4.5 pm. A somewhat narrower subset of both fields was sub- 
sequently covered by Spitzcr-C A N DELS (S-CANDELS; M. 
Ashby et al., in preparation), to even fainter levels; reaching 
~ 0.5 mag deeper than SEDS in both of the warm IRAC band- 
passes. The HUDF09 fields were observed by Spitzer pro- 
gram 70145 (the IRAC Ultra-Deep Field Labbe et al. 2013), 
reaching 120, 50 and 100 hr in the HUDF Main, PARI and 
PAR2 fields, respectively. Finally, program 70204 (PI Fazio) 
observed a region in the ERS field to 100 hr depth. The 
present work is based on mosaics constructed by coadding all 
the above data following the procedures described by Ashby 
et al. (2013). The combined data have a depth of >50 hr over 
both CANDELS GOODS fields and >100 hr over the HUDF 
main field. 

As the IRAC PSF is much broader than that of 1 1ST, our 
galaxies may be blended with other nearby sources. We mea- 
sure Spitzer! IRAC 3.6 and 4.5 pm photometry by perform- 
ing PSF-matched photometry on the combined IRAC data, 
which reach at least 26.8 AB mag (5er) at 3.6 pm (Ashby et 
al. in prep.). We utilized the TPHOT software (Merlin et al. 
in prep.), an updated version of TFIT (Laidler et al. 2007), 
to model low-resolution images (IRAC images) by convolv- 
ing the //-band image with empirically derived IRAC PSFs 
and simultaneously fitting all IRAC sources. Specifically, we 
used the light profiles and isophotes in the detection ( J+H ) 
image obtained by Source Extractor, and then convolved them 
to the high resolution image with a transfer kernel to generate 
model images for the low-resolution data. These models were 
then fit to the real low-resolution images, dilating the segmen- 
tation maps of the model images to account for missing flux 
on the edges of galaxies (Galametz et al. 2013). The fluxes 
of sources are determined by the model which best-represents 
the real data. As the PSF FWHM of the high-resolution im- 
age (//-band) is negligible (~0. 19") when compared to those 
of the low-resolution IRAC images (~1.7"), we use the IRAC 
PSFs as transfer kernels. We derive empirical PSFs by stack- 
ing isolated and moderately bright stars in each field. As our 
own WFC3 catalog was used as the input for TPHOT, all of 
our galaxies have measurements in the TPHOT catalogs. We 
visually inspected the positions of each of our high-redshift 
galaxy candidates in the IRAC images to ensure no significant 
contamination from the residuals of nearby bright galaxies. If 
an object was on or near a strong residual, we ignored the 
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Observed X frim) 

Fig. 4. — The SED of the only galaxy in our 50-object sample of bright 
(Mj 500 < — 21) z > 6 galaxies which had a photometric redshift which pre- 
ferred a low-redshift solution after the inclusion of IRAC and F814W pho- 
tometry. The blue curve shows the original high-redshift best-fitting stel- 
lar population model and photometric redshift probability distribution func- 
tion, while the red curve shows the results including IRAC and F814W. This 
galaxy was removed from our sample, as the IRAC photometry is consistent 
with the stellar emission peak at z ~ 1 . The inferred contamination rate of 
2 % (one out of 50 galaxies) is even lower than our estimates for z > 6 in §3.8. 

IRAC photometry in the subsequent analysis. This was the 
case for 18/164 galaxies at z = 4, 23/85 at z = 5, 3/29 at z = 6, 
6/18 at z = 7 and 1/3 at z = 8. With these contaminated fluxes 
removed, we found that all remaining Misoo < -21 galaxies 
at z= 4-8 had 3.6 pm detections of at least 3cr significance, 
with a magnitude range at z > 6 of 22.7 < m3. 6 < 25.8. The 
full description of our TPHOT IRAC photometry catalog will 
be presented in M. Song et al. (in prep). 

We reran EAZY for this subsample of bright galaxies, in- 
cluding the Spitzer/IRAC fluxes, as well as photometry from 
the ACS F814W filter, which was not included in the orig- 
inal photometric redshift calculation (see §2.1). We exam- 
ined these updated photometric redshift results, searching for 
galaxies in ourz = 4 and 5 samples withz„ ew < 2.5, and in our 
z = 6, 7 and 8 samples with z new < 4. We found 14 out of 164 
galaxies in our z = 4 sample and 14 out of 85 galaxies in our 
z= 5 sample with z new < 2.5. We found one galaxy out of 29 at 
z = 6 that appears to be better fit with a low-redshift solution, 
of z new = 0.9, while zero galaxies in ourz= 7 or 8 samples had 
preferred low-redshift solutions with the inclusion of IRAC 
photometry. 

We examined each of the 29 z = 4-6 galaxies that had pre- 
ferred low-redshift solutions with the inclusion of the IRAC 
photometry to examine whether the updated redshift yielded a 
spectral energy distribution consistent with the observed pho- 
tometry. Out of the 28 z = 4 or 5 galaxies with preferred low- 
redshift solutions, 23 had photometry consistent with a true 
low-redshift galaxy. Four galaxies, however, had photometry 
which appeared to be consistent with a high redshift galaxy 
with a strong emission line (Ha or [O III]) in one IRAC band. 
Systems with lines such as these (i.e., EW[om] >500 A) are 
rare locally, but appear to be more common at high-redshift 
(e.g., van der Wei et al. 2011; Finkelstein et al. 2013; Smit 
et al. 2014). Although typical emission lines strengths are 
now included in the EAZY templates, these do not account 
for extreme emission lines, thus it is not surprising that EAZY 
does not return a high-redshift solution. We elect to keep these 
four galaxies in our sample, noting that the lack of strong 
rest-frame optical lines in the EAZY templates does not af- 
fect our initial sample selection, which does not make use of 


TABLE 2 

Summary of Full Galaxy Sample 


Redshift 

N fl „ 

N m<-21 

Ve/y (NT 1 500 — 22) 
(10 5 Mpc 3 ) 

V e// (M1500 — 19) 
(10 s Mpc 3 ) 

4 (3.5 -4.5) 

4156 

150 

12.2 

4.11 

5 (4.5 -5.5) 

2204 

77 

8.98 

3.36 

6 (5.5 -6.5) 

706 

28 

7.93 

2.50 

7 (6.5 -7.5) 

300 

18 

6.99 

0.30 

8 (7.5 -8.5) 

80 

3 

5.88 

0.16 


Note. — The total number of sources in our final galaxy sample, 
after all contaminants were removed. The final two columns give the 
total effective volume at each redshift for two different values of the 
UV absolute magnitude. 

the IRAC photometry. A fifth galaxy (z5_GNW_13415) has a 
high-quality published spectroscopic redshift ofz= 5.45, thus 
we also keep it in our sample. Perhaps unsurprisingly, the 
low-redshift fit for this galaxy had a very poor quality-of-fit, 
with x 2 = 127, implying that the EAZY templates are a poor 
match for this galaxy. The remaining 23 galaxies at z = 4 and 
5 were removed from our sample. The sole z > 6 galaxy with 
a preferred low-redshift solution with the inclusion of IRAC 
photometry, z6_GSW_3089, is shown in Figure 4. The red 
HST colors imply a much brighter IRAC flux than is seen. 

A solution at z= 0.93 yields a somewhat (though not over- 
whelmingly) better fit, as the peak of stellar emission at that 
redshift better matches the observed IRAC fluxes. We have 
thus removed this galaxy from our sample. We note that this 
galaxy has a spectroscopic redshift of z = 5.59 from the ob- 
servations of Vanzella et al. (2009). However, this object has 
a spectroscopic quality flag of “C”, which indicates that the 
spectroscopic redshift is unreliable. This combined with the 
~ 4<j detection in the F 606 band leaves us confident that a low 
redshift solution is more likely, and removing this object from 
our sample is a conservative approach. 

With the removal of these likely contaminants, we retain a 
total sample of 150, 77, 28, 18 and 3 Mi 5 qo < -21 galaxies at 
z= 4, 5, 6, 7 and 8, respectively. The fraction of contaminants 
at z > 6 (one out of 50 z = 6, 7, and 8 galaxies, or 2%) is 
consistent with (albeit somewhat less than) the expected low 
value calculated in §3.8 below. 

Our final galaxy sample is summarized in Table 2, and a 
catalog of all galaxies in our sample is provided in Table 3. In 
Figure 1 we show the absolute UV magnitude distribution of 
our final samples, highlighting that we cover a dynamic range 
of five magnitudes. In particular, the CANDELS data are cru- 
cial, as galaxies from these data dominate the total number 
of galaxies in our sample, and ^double the luminosity dy- 
namic range which we can probe. This is highlighted in the 
left pane of Figure 2, which shows that galaxies discovered in 
the CANDELS GOODS fields dominate the total number at 
all redshifts in our sample. 

3.7. Stellar Population Modeling 

To derive the rest- frame absolute magnitude at 1500 A 
(A/ 1500 ), as well as the UV spectral slope f3 (fx oc X 13 Calzetti 
et al. 1994), we fit spectral energy distributions (SEDs) from 
synthetic stellar population models to the observed HST pho- 
tometry of our high-redshift candidate galaxies. The tech- 
nique used here is similar to our previous works (Finkelstein 
et al. 2010, 2012b, a, 2013). We used the updated (2007) 
stellar population synthesis models of Bruzual & Chariot 
(2003) to generate a grid of spectra, varying the stellar pop- 
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TABLE 3 

Catalog of Candidate Galaxies at 3.5 <z < 8.5 


Catalog ID 

IAU Designation 

RA 

Dec 

Zphot 

Ml 500 



(J2000) 

(J2000) 


(AB mag) 


z4_GSD_27037 
z4_ERS_3675 
z4_GND_29830 
z4 GND 30689 
z5_GSD 8969 
z5 GND_31173 
z5 MAIN 327 1 
z5_PAR2_3 762 
z6 GND16819 
z6_GNW_16070 
z6_MAIN_29 1 6 
z6_MACS04 1 6PAR145 
z7_GSD_12285 
z7 ERS 6730 
z7_GND_16759 
z7_A2744PAR_4276 
z8 GSD 16150 
z8_MAIN_5173 
z8_GND_32082 
z8 GND 8052 


HRG14 J033240. 8-275003. 1 
HRG14 J033235. 0-2741 17.5 
HRG14 J123718. 1+621309.7 
HRG14 J 12372 1.4+62 1259.2 
HRG14 J0332 1 6.2-274641 .6 
HRG14 J 12373 1.0+62 1254.2 
HRG14 J033243. 5-27471 1.4 
HRG14 J033304. 8-275234. 7 
HRG14 J123718. 8+621522. 7 
HRG14 J123549. 0+62 1224.8 
HRG14 J033244. 8-274656. 8 
HRG14 J04 1 632.2-240533.3 
HRG14 J033206.7-274715.8 
HRG14 J0332 16. 0-2741 59.2 
HRG14 J1236 19.2+62 1523. 2 
HRG14 J00 1357. 5-302358. 3 
HRG14 J0332 13.9-274757.7 
HRG14 J033241. 5-27475 1.0 
HRG14 J 123727.4+62 1244.4 
HRG14 J123704. 8+621718. 8 


53.169922 

53.145882 

189.325211 

189.339355 

53.067379 

189.379272 

53.181351 

53.270004 

189.328232 

188.954025 

53.186806 

64.134117 

53.028114 

53.066677 

189.079834 

3.489512 

53.057983 

53.172874 

189.364258 

189.270020 


-27.834183 3.54(3 

-27.688189 4.08(3 

62.219368 4.01 (3 

62.216450 3.66(3 

-27.778219 5.00(4 

62.215046 4.85 (4 

-27.786510 5.50(4 

-27.876295 4.47 (3 

62.256317 5.55 (5 

62.206898 5.88 (5 

-27.782433 6.42 (5 

-24.092587 5.91 (5 

-27.787714 7.30(6 

-27.699766 6.74 (5 

62.256454 6.69 (6 

-30.399530 6.51 (6 

-27.799349 7.91 (6 

-27.797487 8.11 (6 

62.212334 7.64(7 

62.288559 8.10 (7 


.45 

to 

3.63) 

-20.45 

.79 

to 

4.28) 

-20.39 

.85 

to 

4.17) 

-19.68 

.59 

to 

3.75) 

-21.09 

.87 

to 

5.14) 

-20.62 

.37 

to 

5.09) 

-19.59 

.58 

to 

5.67) 

-16.95 

.71 

to 

4.70) 

-18.83 

.41 

to 

5.65) 

-21.56 

.64 

to 

6.02) 

-20.83 

.79 

to 

6.76) 

-18.39 

.08 

to 

6.23) 

-18.49 

.44 

to 

7.89) 

-19.57 

.64 

to 

6.87) 

-20.31 

.33 

to 

6.89) 

-20.89 

.26 

to 

6.78) 

-19.37 

.21 

to 

8.54) 

-20.14 

.29 

to 

8.64) 

-17.67 

.02 

to 

8.16) 

-20.27 

.04 

to 

8.41) 

-20.68 


(-20.57 to -20.42) 
(-20.47 to -20.15) 
(-19.81 to -19.54) 
(-21.16 to -21.03) 
(-20.68 to -20.51) 
(-19.77 to -19.43) 
(-17.05 to -16.78) 
(-18.97 to -18.57) 
(-21.62 to -21.47) 
(-20.88 to -20.64) 
(-18.55 to -18.19) 
(-18.64 to -18.16) 
(-19.79 to -19.31) 
(-20.31 to -19.96) 
(-20.98 to -20.76) 
(-19.47 to -19.22) 
(-20.38 to -19.91) 
(-17.98 to -17.41) 
(-20.33 to -20.02) 
(-20.84 to -20.44) 


Note. — A catalog of our 7446 z = 4-8 galaxy candidates, with their derived properties. We include the IAU designation 
for continuity with previous and future works, with a designation prefix HRG14 denoting "High Redshift Galaxy 2014”. 
The values in parentheses represent the 68% confidence range on the derived parameters. Here, we show 20 representative 
galaxies, four from each redshift bin. This table is available in its entirety in a machine-readable form in the online journal. A 
portion is shown here for guidance regarding its form and content. 


illation metallicity, age, and star-formation history. Metallic- 
ities spanned 0.02 -lx Solar, and ages spanned 1 Myr to 
the age of the universe at a given redshift. We allowed several 
different types of star-formation histories (SFHs), including a 
single burst, continuous, as well as both exponentially decay- 
ing and rising (so-called "tau” and "inverted-tau” models). To 
these spectra, we added dust attenuation using the starburst 
attenuation curve of Calzetti et al. (2000), with a range of 0 < 
E(B-V) < 0.8 (0 < A v < 3.2 mag). We also included nebular 
emission lines using the prescription of Salmon et al. (2014), 
which uses the line ratios from Inoue (2011), based on the 
number of ionizing photons from a given model, and assum- 
ing the ionizing photon escape fraction is « zero. We then 
redshifted these models to 0 < z < 11 and added intergalac- 
tic medium (IGM) attenuation (Madau 1995). These model 
spectra were integrated through our HST filter bandpasses to 
derive synthetic photometry for comparison with our obser- 
vations. For each model, we computed the value o f M ] 500 by 
fitting a 100 A-wide synthetic top-hat filter to the spectrum 
centered at rest- frame 1500 A. Likewise, for each model we 
measured the value of / 3 by fitting a power law to each model 
spectrum using the wavelength windows specified by Calzetti 
et al. (1994), similar to Finkelstein et al. (2012b). 

The best-fit model was found via y 2 minimization, includ- 
ing an extra systematic error term of 5% of the object flux 
for each band to account for such items as residual uncer- 
tainties in the zeropoint correction and PSF-matching pro- 
cess. The stellar mass was computed as the normalization 
between the best-fit model (which was normalized to 1 M 0 ) 
and the observed fluxes, weighted by the signal-to-noise in 
each band. These best-fit values of AT 1500 are used in our lumi- 
nosity function analysis below, while j3 is used to correct for 
incompleteness in a color-dependent fashion. The uncertain- 
ties in the best-fit parameters were derived via Monte Carlo 
simulations, perturbing the observed flux of each object by a 
number drawn from a Gaussian distribution with a standard 


deviation equal to the flux uncertainty in a given filter. For 
each galaxy, 10 3 Monte Carlo simulations were run, provid- 
ing a distribution of 10 3 values for each physical parameter. 
The 68% confidence range for each parameter was calculated 
as the range of the central 68% of results from these simula- 
tions. In these simulations, the best-fit redshift was allowed 
to vary following the redshift probability distribution func- 
tion, thus folding in the uncertainty in redshift into the un- 
certainty in the physical parameters (most notably, the stellar 
mass andMisoo; Finkelstein et al. 2012a). During this process, 
we only allowed the redshift to vary within Az = ± 1 of the 
best-fit photometric redshift. This excludes any low-redshift 
solution from biasing the uncertainties on a given parameter. 
The amount of the integrated P(z) at z > 3 excluded via this 
step was typically < 10% (at z = 6). 

3.8. Contamination 

A key issue in any study of high redshift galaxies is the 
risk of sample contamination, either by spurious sources or 
by lower-redshift interlopers. The gold standard for eliminat- 
ing contamination is to obtain spectroscopic redshifts. This is 
clearly unfeasible for all galaxies in our sample (until the next 
generation of space and ground-based telescopes), but there 
is significant archival spectroscopic data. As discussed in 
§3.2, we find excellent agreement between available spectro- 
scopic redshifts and our photometric redshifts, with (J z /(\+z) = 
0.031. In particular, the four bright 4 galaxies (24. 9 < A 25 < 
25.7) with confirmed z spec >6.5 have quite excellent agree- 
ment: z7_GNW_24443 with z p h ot = 6.66 and z spec = 6.573 
(Rhoads et al. 2013), z7_GSD_21172 with z p / wt = 6.73 and 
Zspec = 6.70 (Hathi et al. 2008), z7_GNW_4703 with Z phot 

4 Object z7_MAIN_2771 has a tentative spectroscopic redshift of z spec = 
7.62, based on a 4a possible Ly a emission line from Schenker et al. (2014). 
This object is quite faint, with J 125 = 28.8, resulting in a somewhat large 95% 
photometric redshift confidence range of 6.02-7.74, though consistent with 
the tentative spectroscopic redshift. 
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z=6 z=7 z=8 



Z Z Z 


Fig. 5. — Redshift PDFs for galaxies in each of our three redshift samples, stacked in bins AMj/f = 1. The legends give the number of galaxies in each stack, 
as well as the fraction of the redshift PDF at z > 4 (denoted as P). Even in the worst case (which is for faint-galaxies at z = 6) <16.3% of the sample could 
possibly be at lower redshift. 


7.19 and z spec = 7.213 (Ono et al. 2012), and z7_GND_42912 
with z p h 0 t = 7.45 and z spec = 7.51 (Finkelstein et al. 2013). 5 
The brightest source in our z= 6 sample, z6_GSW_12831 
with A /1500 = -22. 1 and z p ] lot = 5.77, is confirmed with z spec = 
5.79 (Bunker et al. 2003). This galaxy has a 3 a detection in 
the F 606 -band, which could have resulted in its exclusion from 
a typical LBG color-color selection sample (e.g., Bouwens 
et al. 2007), though the observed flux can be explained by 
non-ionizing UV photons transmitted through the Lya forest. 

In general, spectroscopic followup of sources selected on 
the basis of their Lyman breaks (either color-color selection, 
or photometric redshift selection) finds a very small contam- 
ination by low-redshift sources (e.g., Pentericci et al. 2011). 
However, given the apparent difficulty in detecting Lya emis- 
sion at z > 6.5 (e.g., Pentericci et al. 2011; Ono et al. 2012; 
Finkelstein et al. 2013), the true effect of contamination at 
these higher redshifts is not empirically known. In this sub- 
section, we will attempt to estimate our contamination frac- 
tion by other means. 

3 . 8 . 1 . Properties of the Image Noise 

Two key components of our selection processes should 
eliminate contamination by spurious sources in our sample. 
First, we restricted our sample to galaxies detected at >3.5<r 
in two imaging bands: J 125 and Formally, requiring a 

3.5 a detection in a single band should yield only a 0.05% con- 
tamination by noise. However, the wings of the noise distribu- 
tion are highly non-Gaussian. We examined this by measuring 
the fluxes at 2 x 10 5 random positions in the J 125 and//i 6 o im- 
ages in the GOODS-S Deep field (see Schmidt et al. (2014) 
for a similar analysis). To avoid biasing from real objects 
in the image, we only considered negative fluctuations (e.g., 
Dickinson et al. 2004), where the contamination percentage 
was computed as the ratio of the number of apertures with a 
flux < -1x3. 5cr to the total number of apertures with neg- 
ative fluxes. We found that in each of these bands individu- 

5 This source was originally called z8_GND_5296 in our previous cata- 
log (Finkelstein et al. 2013). Our new catalog uses an updated version of the 
CANDELS GOODS-N data, thus the catalog numbering is different. Addi- 
tionally, the slightly updated photometry pushes the photometric redshift of 
this galaxy slightly below z = 7.5, placing it in the z = 7 sample. 


ally, the fraction of positions measuring at >3.5er was ~1.4%, 
much higher than predicted based on an assumption of Gaus- 
sian noise. If we instead look at the number of 3.5<r fluctu- 
ations in both the J \ 25 and / f <-,o images at the same location, 
we find that this contamination drops to nearly zero, at 0.05%. 
Thus, we conclude that as we require significant detections in 
two bands, the contamination in our sample by noise is negli- 
gible. Spurious sources other than noise spikes are eliminated 
by our detailed visual inspection of each source, described in 
§3.3. 

3 . 8 . 2 . Estimates from Stacked Redshift PDFs 

Contamination by low-redshift interlopers is a more com- 
plicated issue. While extreme emission line galaxies at lower- 
redshift could theoretically be an issue (Atek et al. 2011), our 
requirement of detections in two bands (as well as the fre- 
quent detections in more than two bands for all but the highest 
redshift objects in the z= 8 sample) makes a significant con- 
tamination by these sources unlikely. The most likely possible 
contaminants are faint red galaxies at z < 2 (e.g., Dickinson 
et al. 2000). These galaxies can be too faint to be detected in 
our optical imaging, but their red SEDs yield detections in the 
WFC3/IR bands. Although faint sources that are very red will 
have a disfavored high-redshift solution with our current pho- 
tometric selection, we have information on their likelihoods 
encoded in our redshift PDFs. Figure 5 shows the redshift 
PDFs of galaxies in each of our three highest redshift samples, 
stacked in magnitude bins of AM = 1 mag. At all redshifts 
and all magnitudes, >85% of the redshift PDF is at z > 4, 
implying that there is not significant contamination by lower- 
redshift galaxies. The position of the secondary redshift peak 
is consistent with the redshift were our detected spectral break 
due to a 4000 A break rather than the Lyman break (at z = 6 , 
7 and 8 , this gives z secon d aiy = 1.1, 1,4 and 1.7). At z= 8 , the 
possible contamination is <10.5%, primarily due to the fact 
that at z = 8 , a galaxy will have to be undetected in most of 
the filters we consider here. The worst case is for faint galax- 
ies at z = 6 , as z = 6 galaxies are typically detected in all but 
two filters, though even here, the indicated contamination is 
<15%. 
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Fig. 6. — Top: Filter transmission curves for the filter set used in this study (T 098 , which was used in the GOODS-S ERS field only, is not shown). The vertical 
lines denote the relative position of the Lya break (rest 1216 A) in a given filter for galaxies at the center of our three highest redshift bins. Bottom) Image 
stacks of galaxies in our three highest redshift galaxy samples. If our sample had a significant fraction of lower-redshift interlopers, significant flux would be 
seen blueward of the break (e.g., .#435 at z = 6, z'775 and blueward at z = 7 , and 7 g 14 and blueward at z = 8). This is not observed at any redshift, thus we conclude 
that our sample is free of significant contamination. 


3.8.3. Stacking Imaging 

The limits from the previous subsection are likely upper 
limits on the contamination fraction. When fitting photomet- 
ric redshifts, to rule out all low-redshift solutions, the Lyman 
break needs to be detected at high significance, which is the 
case for only the brightest galaxies (e.g., at z = 6, the brightest 
bin has a contamination of <2%). Additionally, these results 
are dependent on the templates used, which by definition do 
not account for unknown galaxy populations. We therefore 
consider two empirical tests of contamination. The first is to 
stack all galaxies in a sample, in order to search for detections 
below the Lyman break. The results from this test for z = 6, 
7 and 8 are shown in Figure 6. As expected for galaxies at 
the expected redshifts, there is no visible signal in the B 435- 
band at z = 6, «775-band and blueward at z = 7, and /si4- band 
and blueward at z = 8. This confirms our previous conclusion 
that there is not a significant contamination by lower-redshift 
sources in our sample. 

3.8.4. Estimates from Dimmed Real Sources 


a given field with a similar magnitude as the dimmed source. 
We then added scatter to the dimmed fluxes, perturbing them 
by a random amount drawn from a Gaussian distribution with 
a standard deviation equal to the flux uncertainties of the ob- 
ject. We included two realizations of the HUDF field to in- 
crease the number of dimmed objects. 

The total number of sources in our artificially dimmed cat- 
alog was 4066 in the Deep fields, and 1254 in the FlUDF 
field. We measured photometric redshifts of these sources 
with EAZY in an identical manner as on our real catalogs, 
and then we applied our sample selection to this dimmed cat- 
alog. In the Deep fields, we found a total number of 149, 134, 
54, 23 and 8 dimmed objects satisfied our z = 4, 5, 6, 7 and 
8 selection criteria. Investigating the original (not dimmed 
or perturbed) colors of these sources, we found that they are 
unsurprisingly red, with the bulk of sources having V — H > 
2 mag. It is therefore this parent population of red sources 
which are responsible for the majority of the possible con- 
tamination. The contamination fraction in our high-redshift 
sample was then defined as 


As a final test, we estimated the contamination by artifi- 
cially dimming real lower redshift sources in our catalog, to 
see if the increased photometric scatter allows them to be se- 
lected as high redshift candidates. This empirical test is use- 
ful as it does not rely on known spectral templates to derive 
the contamination, though it does assume that the fainter ob- 
jects which could potentially contaminate our sample have 
similar SEDs to the bright objects which we dim. We per- 
formed this exercise twice, once using the combined catalog 
of the GOODS-S and GOODS-N Deep fields, and once in 
the HUDF main field, to probe fainter magnitudes. In the 
GOODS Deep fields, we selected all real sources with 21 
< //160 < 24 and z p h ot < 3, and reduced their observed fluxes 
by a factor of 20. The same was done for sources drawn from 
the HUDF Main field, here extending the magnitude range to 
be 21 < //160 < 26. The limits on these magnitudes were cho- 
sen to exclude any real high-redshift sources. We replaced the 
true flux uncertainties of these objects with flux uncertainties 
from a randomly drawn real source from the full catalog from 


T = 


N d immed,sel ect 
N dimmed, red 


* Bl total, red 


N z 


( 1 ) 


where, N dimmed, select was the number of dimmed sources satis- 
fying our high-redshift sample selection, Ndimmed,red was the 
total number of sources in the dimmed catalog with origi- 
nal colors of V — H < 2, N totalled was the number of sources 
in the full object catalog with 25 < H < 27, z p i, ot < 3, and 
V—H < 2, and N z was the number of true galaxy candidates in 
a given redshift bin. For example, at z= 6, where we found 54 
dimmed galaxies satisfied our selection criteria (N dimmed,. select 
=54), N dlmmed%red = 1023, N tot aired = 695, and N-_ = 322, giving 
an estimated contamination fraction T = 11.4%. Therefore, 
for sources with 25 < H < 27, we found an estimated con- 
tamination fraction of T = 4.5%, 8.1%, 11.4%, 11.1% and 
16.0% at z = 4, 5, 6, 7 and 8, respectively. We performed the 
same exercise in the HUDF, here for fainter sources with 26 
< H < 29, finding 30, 21, 8, 8 and 0 sources satisfied our z = 
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4, 5, 6, 7 selection criteria, giving a contamination fraction of 
9.1%, 11.6%, 6.2%, 14.7% and <4.9% atz = 4, 5, 6, 7 and 8. 

Broadly speaking, we estimate relatively small contamina- 
tion fractions of ~5-15%, in-line with the estimates above 
from the stacked P(z) curves. As the bulk of contaminants 
appear to be red galaxies, it is interesting to compare to the 
space density of these potentially contaminating sources. This 
was recently estimated by Casey et al. (2014b), who find that 
dusty star-forming galaxies at z < 5 will contaminate z > 5 
galaxies samples at a rate of <1%. This is much less than our 
contamination estimates, thus we may have overestimated the 
contamination rate, though it may not be inconsistent once 
photometric scatter is applied to faint, red galaxies, making it 
easier for them to scatter into our sample. In any case, the ex- 
pected contamination rate is quite small, therefore we do not 
reduce our observed number densities for the expected mini- 
mal contamination. 

4. COMPLETENESS SIMULATIONS 

We performed an extensive set of simulations to estimate 
the effective volume for each source in our sample, account- 
ing for both image incompleteness and selection effects. We 
inserted mock galaxies into the imaging data, repeating the 
same analysis for source detection, photometry, photometric 
redshift measurement, and sample selection as was done on 
the real data. We then compared the fraction of recovered and 
selected mock sources to the total number of input sources in 
a given bin of absolute magnitude and redshift to determine 
our completeness in that bin. 

While the effective volumes are typically computed as a 
function of magnitude and redshift, other key factors in these 
simulations are the choices of galaxy size and color. At a con- 
stant magnitude, a very extended or very red galaxy may not 
make it into our sample, thus the effective volume depends 
not only on magnitude and redshift, but also on the size and 
rest-frame UV color. To see what effect this has, we have 
computed our completeness as a function of four properties: 
redshift, absolute magnitude, half-light radius, rest-UV color, 
where we have parametrized the latter via the UV spectral 
slope f3 (Calzetti et al. 1994). A large number of simulated 
objects are needed to fill out this four-dimensional space; our 
completed simulations recovered ~5.4 million out of 7 mil- 
lion objects input across all of our fields (where the recovered 
objects were detected in our photometry catalogs; this number 
does not account for the photometric redshift selection, which 
we discuss below). 

Our simulations were run separately on each of our 10 sub- 
fields defined in Table 1. To ensure that the mock galax- 
ies did not affect the background estimation, a small num- 
ber of galaxies were added during each simulation. To opti- 
mize the simulation runtime, the mock galaxies were added to 
cutouts from the full images. In the GOODS sub-fields (i.e., 
CANDELS Deep and Wide, and the ERS), 200 mock galax- 
ies were added to a 2000x2000 pixel (2' x 2') region of the 
images, while for the single-pointing FlUDF and HFF fields, 
100 galaxies were added to a 1000x1000 pixel region. As 
the depth across our imaging data can vary, during each sim- 
ulation the position of the cutout varied, such that when we 
combine all of our simulations, we average over any differ- 
ences in the depth across a given field. 

To determine the colors of the mock galaxies, we created 
distributions in redshift, dust attenuation (parameterized by 
E[B-V]), stellar population age and stellar metallicity. The 
redshift was defined to be flat across 3 < z < 9, such that we 


simulate objects well above and below the redshift ranges of 
interest. The dust attenuation E(B-V) was defined to have a 
Gaussian distribution with a mean of 0. 1 and a a = 0. 15 (with 
a minimum of zero). The age was defined as a log-normal dis- 
tribution, with a peak near 10 Myr, and a tail extending out to 
the age of the universe at a given redshift. The metallicity dis- 
tribution was also log-normal, with a peak of Z = 0.2 Zq, and 
a tail towards higher values. The exact values of these dis- 
tributions are not crucial given our methodology (as opposed 
to a multivariate analysis, where the distributions are very im- 
portant), as they combine to create a distribution of rest-frame 
UV slope f3. We crafted these distributions to provide a dis- 
tribution of (3 encompassing the expected values for our real 
objects. We then used the updated (2007) stellar population 
models of Bruzual & Chariot (2003) to calculate the colors of 
a stellar population given the distributions above. To convert 
these colors into magnitudes, we assumed a distribution of H- 
band magnitudes designed to have many faint ( H > 26) galax- 
ies (which is where we expect to become incomplete), and 
relatively few at bright magnitudes. To ensure enough bright 
galaxies to calculate a robust incompleteness, every 10th sim- 
ulation used a flat distribution of //-band magnitudes of 22 
< H < 25. These H - band magnitudes were combined with 
the mock galaxy colors to generate magnitudes in each filter 
for a given field. 

To generate the galaxy images themselves, we used the 
GALFIT software (Peng et al. 2002). We assumed a log- 
normal distribution of half-light radii with a peak at 1-pixel, 
and a high tail towards larger radii, giving an interquartile 
range of half-light radii of 1.4-4. 9 pixels. This corresponds 
to ~0.4-1.6 kpc, spanning the range of the majority of re- 
solved galaxies at z > 4 (e.g., Oesch et al. 2010a; Grazian 
et al. 2012; Ono et al. 2013; Curtis-Lake et al. 2014). GAL- 
FIT also requires a Sersic index («), axis ratio and position 
angle; the Sersic index was assumed to be a log-normal dis- 
tribution at 1 < n < 4, with the majority of the mock galaxies 
having disk-like morphologies (« < 2); the axial ratio was also 
log-normal, with a peak at 0.8, and a tail toward lower values, 
and the position angle was a uniformly distributed random 
value between 0 and 360 degrees. GALFIT was then used to 
generate a 101 x 101 pixel (6" x 6") stamp for a given mock 
galaxy, which was then added to the image at a random loca- 
tion. Because our data are PSF-matched to the H- band, we 
had GALFIT convolve the mock galaxy images with the H- 
band PSF prior to adding them to the data for all filters. 

Once the set of mock galaxies for a given simulation were 
added to the data, photometric catalogs were generated us- 
ing Source Extractor in the exact same manner as was done 
on the data. This requires making the same set of detection 
images, which we did here using weighted combinations of 
the simulated images. These catalogs were read in and com- 
bined, again in the same methodology as with the data, in- 
cluding aperture corrections (the exception here is that a cor- 
rection for Galactic extinction was not applied in the simu- 
lation, as the simulated objects did not have Galactic extinc- 
tion included). The photometric catalog was then compared 
to the input catalog to generate the list of recovered objects 
(i.e., mock galaxies which were recovered by Source Extrac- 
tor); an object was regarded as being recovered if it had a 
positional match within 0.2" of one of the input mock galax- 
ies. The recovered object catalogs were processed through 
EAZY to generate photometric redshifts, and then run through 
our SED-fitting routine to measure absolute UV magnitudes 
(A/i 5 oo), stellar masses and UV spectral slopes. These simu- 
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Fig. 7. — The results of our completeness simulations, showing the probability that a given simulated source was recovered as a function of its input redshift. 
The solid lines denote sources with M 1500 = — 22, while the dashed lines denote M 1500 = _ 19. These lines assume a half-light radius of r /, = 0. 18 /; and (3 = —2.0. 
The background histogram shows the distribution of best-fit photometric redshifts for the real galaxies in each redshift subsample. Although our selection criteria 
combined with the wavelengths probed by our filter-set results in a completeness that peaks at close to z = 4, 5, 6, 7 or 8, the evolving luminosity function as well 
as our sensitivity to bright galaxies results in our samples having mean redshifts slightly lower than the bin center, particularly in the higher redshift samples. 


lations were then repeated until a large sample of recovered 
galaxies was available, which were then compiled in a single 
database per field. The completeness was defined as the num- 
ber of galaxies recovered versus the number of input galaxies, 
as a function of input absolute magnitude, redshift, half-light 
radius and UV spectral slope {3. Figure 7 shows the results 
from our simulations. 

In our original simulations the recovered redshift was typi- 
cally ~0.2 lower than the input redshift, independent of mag- 
nitude. This is likely not a fault in our photometric redshift 
estimates, as Figure 2 shows that these agree well with ex- 
isting spectroscopic redshifts for real galaxies. Rather, it is 
likely a mismatch between our simulated SEDs and those of 
the templates used in EAZY. Upon further investigation, we 
found that the cause of this offset was Lyo emission in the 
mock galaxies. While Lya photons were attenuated by dust 
in the same manner as adjacent UV photons, we did not in- 
clude any additional Lya attenuation for, e.g., geometric or 
kinematic effects. This led to very high Lya escape frac- 
tions, which were not matched in the templates. This high 
Ly a emission reduced the amplitude of the Lyman break, re- 
sulting in a (slightly) lower photometric redshift. After re- 
ducing the amount of Lya flux to 25% of the intrinsic value, 
our photometric redshifts matched the input redshifts. Rather 
than rerun all of our completeness simulations, we elected to 
simply reduce the input redshift by 0.2 when interpreting our 
simulations, which corrects for this effect (this changes the 
distance modulus by <0.1 mag). The exception was the sim- 
ulations for the HFF parallel fields, which were run after this 
effect was noticed. In those fields, the input models had their 
Ly a flux reduced to 25% of the intrinsic value, and no change 
to the model redshift was needed. 

It is important to examine whether the choice of comput- 
ing the completeness as a function of input properties affects 
our result. As mentioned in §2.3, we used the results from 
these simulations to correct for offsets in the recovered versus 


input magnitudes (i.e., to be sure the fluxes we use represent 
the total flux). Additionally, we examined whether there ex- 
ist biases in the half-light radius or [3 measurements from the 
simulations. Recovered objects were typically measured to 
have a half-light radius ^0.03" (0.5 pixels) smaller than the 
input value, and were measured to be slightly redder (A/3 < 
0.1). However, these corrections make effectively no change 
to the effective volumes derived from the simulations, and so 
were not applied. 

In each redshift bin, the effective volume for galaxies in a 
given field was then calculated via 

fdV 

V e ff (M 1500 ,r h ,[3ym / — P(M 1500 ,z,r h , (3) dz (2) 

where dV /dz is the comoving volume element, and 
P(Mi 5 oo,z,rh,pi) is the result from our completeness simu- 
lations. The integral was done over A z = 1, centered on 
the center of each redshift bin. In each field, we used a 
weighted mean of this three-dimensional effective volume 
V e ff(M 1 500 , Ai , 0) to calculate V eff (M l500 ), where the weight- 
ing is based on the number of real objects in a magnitude bin 
with a given value of r/, and (3, as 


V e ff(M 15 oo) 


E V e fj{M- 1 500 , fh , /?) N(M 1500 , r / t , (3) 

rh 0 

rh 0 


(3) 


This assumes that the completeness corrections estimated us- 
ing our observed size and color distributions are similar to 
what we obtained if we could measure the true sizes and col- 
ors, motivated by our measurement of minimal size and color 
biases when comparing the input to recovered values. 

This weighted volume is the most representative of the true 
volume we are sensitive to, as we explicitly account for the 
incompleteness as a function of size and color. Figure 9 high- 
lights the dependence of the effective volume on these quan- 
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Fig. 8. — The effective volume per unit area of our survey for high-redshift galaxies in each of our redshift bins. Here we divide out the area of a given 
sub-field such that one can easily compare the completeness as a function of magnitude of the various fields. The solid lines give way to dashed lines when the 
volume per unit area falls below 50% of the maximum value. For the luminosity function, we only consider magnitude bins in each field brighter than these 50% 
completeness points, to avoid having data dominated by incompleteness corrections. We note that the ERS has a different 7-band filter (7o98), which gives a 
better spectral resolution around 1 fim, likely responsible for the increased selection efficiency at z = 8 in that field. 



Fig. 9. — Left) The effective volume as a function of UV absolute magni- 
tude for galaxies in our z = 6 sample in the GOODS-N Deep field. The red 
line shows the mean effective volume for this field, weighted by the number 
of galaxies at a given radius and UV color. The black lines show how the 
effective volume changes as a function of effective radius (r e ) for a fixed UV 
color Q 3 = -2). Our weighted mean volume is similar to the effective volume 
assuming r e = 2.5 kpc for bright galaxies, and r e = 1.1 kpc for faint galaxies. 
Middle) The dependence of the effective volume on effective radius in two 
magnitude bins. At fainter magnitudes, the effective volume drops steeply 
with increasing size, as the surface brightness drops below detectable levels. 
Right) Same as middle, except here showing the dependence on UV color. 
The dependence on color is much weaker than that on size, as we remain sen- 
sitive to galaxies until (3 becomes redder than -1, which is much redder than 
the colors of observed high-redshift galaxies (e.g., Finkelstein et al. 2012b; 
Bouwens et al. 2013) 

tities for galaxies in our z = 6 sample in the GOODS-N Deep 
field. The effective volume has a strong dependence on the 
surface brightness of galaxies, as the volume drops steeply 
both for larger sizes and fainter magnitudes. The central and 
right panels highlight that while the effective volume (and 
thus sample completeness) is sensitive to both size and color, 
the color has a relatively minor role. We remain sensitive to 
fairly red galaxies {(3 = -1), similar to previous results from 
Bouwens et al. (2012). We note that although the effective 
volume has a strong dependence on size, the relatively small 
sizes of galaxies in our sample yields a volume similar to 
that obtained when assuming a constant small size. Thus, al- 
though our volumes are the most accurate, had we assumed a 


fixed effective radius of, e.g., r e = 1 kpc, our results would not 
change significantly. This is consistent with the conclusions 
of Grazian et al. (2012) who found, accounting for the size- 
luminosity relation when deriving the z = 7 luminosity func- 
tion, similar results as previous studies that neglected the size- 
luminosity relation. Our final effective volumes are shown in 
Figure 8. 

5. THE LUMINOSITY FUNCTION 
5.1. Parametric Approach 

Possessing our final galaxy sample with measured values 
of Mi 5 oo, as well as the effective volumes from our complete- 
ness simulations in the previous section, we can now proceed 
to measure the rest-frame UV luminosity function at z = 4, 
5, 6, 7 and 8. We calculate the luminosity function in two 
ways: a parametric version assuming that the luminosity func- 
tion takes the form of a Schechter (1976) function, and a non- 
parametric step-wise maximum likelihood (SWML) calcula- 
tion. 

The fitting of a Schechter function is well motivated, as it 
successfully matches the observed rest-UV luminosity func- 
tions at lower redshifts (e.g., Reddy & Steidel 2009; Bouwens 
et al. 2006). This function is characterized by a power-law 
at the faint-end with slope a, and an exponential cut-off at 
the bright end, transitioning between the two regimes at the 
characteristic magnitude M*. The parameter p* sets the nor- 
malization of this function. The number density at a given 
magnitude is then given by 

p(M) = 0.4 In (10) (4) 

For the measurement of the luminosity function assuming 
a Schechter functional form, we calculated the likelihood that 
the number of observed galaxies in a given magnitude bin is 
equal to that for an assumed value of the Schechter parameters 
M* and a. Rather than performing a grid-based search, we 
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performed a Markov Chain Monte Carlo (MCMC) search al- 
gorithm, to better span the parameter space, as well as to bet- 
ter characterize the uncertainties on the Schechter parameters. 
We performed this calculation in bins of absolute magnitude 
with AM = 0.5 mag, ranging from -24 < M1500 < -17. At the 
bright end we are in the limit of small numbers, therefore we 
model the probability distribution as a Poissonian distribution 
(e.g., Cash 1979; Ryan et al. 2011), with: 

C\v) = ~2\nC(v) (5) 

C (V?) ^ 'y ^ ^ A/j, 0 bs model) — ^Yj, model — ln(-^,obs • ) (6) 

< j 

where C{p) is the likelihood that the expected number of 
galaxies (A mo d e i) matches that observed (IVobs) for a given 
value of M* and a, and C 2 is the goodness-of-fit statistic. 
The subscripts i and j represent the sub-fields and magnitude 
bins, respectively. The final goodness-of-fit is the sum over 
all fields and magnitudes in a given redshift bin. We use the 
effective volume results for a given redshift, magnitude bin, 
and field to convert from the model number density to the ex- 
pected number, calculating p* as the normalization such that 
the total expected number of galaxies over all magnitude bins 
matches the total number of observed galaxies. 

For each magnitude bin, we performed 10 independent 
MCMC chains utilizing a Metropolis-blastings algorithm, 
each of 10 5 steps, building a distribution of M*, a and tp* 
values for each field. During each step of the chain, the likeli- 
hood of a given model was computed for each of our observed 
fields, and then added together to compute the likelihood for 
the sample as a whole (we also recorded the individual field 
values, see §6.5). Prior to each recorded chain, we performed 
a burn-in run with a number of steps equal to 10% of the num- 
ber in each chain. The starting point for the burn is a brute- 
force x 2 fit °f a grid °f a an d M* values to our data. At the 
end of the burn, the final values of the parameters from the last 
step were then the starting points for each chain. The burn-in 
results were not otherwise recorded. During each step, new 
values of M* and a were chosen from a random Gaussian 
distribution, with the Gaussian width tuned to generate an ap- 
proximate acceptance rate of 23%. During each step p* was 
calculated as the normalization. If the difference between the 
likelihood of the model for the current step exceeds that from 
the previous step by more than a randomly drawn value (= 
2 In (/?); where n is a uniform random number between zero 
and unity), then the current values of the Schechter function 
parameters were recorded. If not, the chain reverted to the 
value from the previous step. 

By running 10 independent chains, we mitigate against be- 
ing trapped by local minima in the parameter space. Our final 
result concatenates these 10 chains together, giving a distribu- 
tion of 1 0 6 values of the Schechter function parameters at each 
redshift. The results were visually inspected to confirm that 
the chains reached convergence. For each Schechter function 
parameter, the best-fit values were taken to be the median of 
the distribution, with the uncertainty being the central 68% of 
the distribution. These results are given in Table 4. For the z = 
8 Schechter function fit, we imposed a top-hat prior forcing 
My v to be fainter than -23. Without this prior, the fit pre- 
ferred a much brighter value of /l/( JV , such that the observed 
data points all lay on the faint-end slope (i.e., a single power 
law). We discuss the implications of this in §6.6. 

Although we computed the volumes down to very faint 


TABLE 4 

Schechter Function Fits to the Luminosity 
Function 


Redshift 

M* 

a 

(Mpc~ 3 ) 

4 

-20.73!°;°:; 

-1 56 +006 
- 0.05 

(14.l!f-g|) xltr 4 

5 

— 20.81!°;!° 

-1 67 +005 
^'-o.oo 

(8.95!|;°j) x 10^ 

6 

-2i n +0 - 25 

- 0.31 

—2 02 +01 ° 
- 0.10 

(1.86!°$) x IO - 4 

7 

-21 03 +0 - 37 
Z1 - - 0.50 

-2 03 +0 - 21 
^ - 0.20 

(1.57!>;«) xlO- 4 

8 

-20 89 +0 - 74 
zu - - 1 .08 

-2 36 +0 - 54 
- 0.40 

(0.72! 2 0 ;g) xliT 4 


Note. — The final values for each parameter 
are the median of the parameter distribution from the 
MCMC analysis. The quoted errors represent the 68% 
confidence range on each parameter. 

magnitudes, should we include these faint galaxies and calcu- 
late the luminosity function down to A/1500 =-17 or fainter, we 
would be highly incomplete (Figure 8). In practice, it is 0111- 
deepest field (the HUDF) which determines how faint we can 
constrain the luminosity function. The HUDF drops below 
50% completeness at magnitudes fainter than M1500 ~ -17.5 
at z = 4, 5 and 6, -18 at z = 7, and -18.5 at z = 8. Thus, in 
our calculation of the luminosity function, we only include a 
given field’s contribution at a given magnitude if it is above 
the 50% completeness limit for that magnitude and redshift. 
Extending the analysis fainter will give results dominated by 
the incompleteness correction. 

As shown in Figure 8, while the volume per unit area for the 
different fields is very tight at z = 4, there is a progressively 
larger scatter apparent when moving towards higher redshift, 
representing a systematic uncertainty in the effective volume 
calculation. One likely culprit is the fact that the volumes de- 
pend on the distribution of the sizes and colors of objects in a 
given field. For fields with few sources (i.e., the smallest fields 
at the highest redshifts), there may be only a single object in a 
given magnitude bin. To mitigate significant variances in the 
effective volume at the bright end, where numbers of sources 
are small, we set the effective volume in a given redshift bin 
and field in bright bins with less than three objects equal to the 
value in the brightest bin with more than three objects (i.e., if 
there are no magnitude bins with more than three objects at 
M < - 21, the effective volumes for all brighter bins are set 
equal to the value at M < - 21). This change has no discern- 
able effect on our luminosity function results as this is well 
above the 90% completeness limit for any of our fields, and 
is thus only done to keep small numbers of galaxies from sig- 
nificantly affecting the volumes. 

Another possible issue is the source density of simulated 
objects. In the smaller fields (HUDF, HUDF parallels, HFF 
parallels) we input sources with twice the surface density as 
in the larger fields to speed up the computing time. As some 
sources (just like real galaxies) will inevitably fall on top of 
real sources, and thus not be recovered, an increased source 
density could result in a (slightly) lower completeness. This is 
just what is observed in these fields, as shown in Figures 7 and 
8. To account for this uncertainty, we measured the spread 
in volume per unit area in each field at M1500 = -21 at each 
redshift, which we found to be ~1.5%, 3.8%, 6.2%, 7.8% 
and 13% at z = 4, 5, 6, 7 and 8, respectively. At each step 
in the MCMC chain, we perturbed the effective volume by 
this amount to account for this systematic uncertainty in our 
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Fig. 10. — The rest-frame UV luminosity functions for our z = 4-8 galaxy samples. The large red circles denote our step-wise maximum likelihood luminosity 
function, while the solid red line denotes our best-fitting Schechter function, with the best-fit values given by the inset text. We do not make use of our data below 
the determined 50% completeness level in each field. As the HUDF is our deepest field, the magnitude of our last data point denotes the 50% completeness limit 
in the HUDF. The dashed line shows the best-fit single power law. We also show several luminosity functions from the literature as indicated in the legends. 


luminosity function results. 

5.2. Non-Parametric Approach 

We have also examined a non-parametric approach to 
studying evolution in the luminosity function. This is par- 
ticularly warranted at very high redshift, where the effects 
responsible for suppressing the bright-end of the luminosity 
function and causing the exponential decline in number den- 


sity (e.g., active galactic nuclei feedback, or dust attenuation) 
may be less relevant. We thus calculated the SWML luminos- 
ity function, which is essentially the number density at a given 
magnitude, free from assumptions about the functional form 
of number density with magnitude. We also calculated the 
SWML luminosity function using an MCMC sampler. In this 
case, as the number densities in the magnitude bins are not 
linked by an overarching function, we calculate the number 
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TABLE 5 

Rest-Frame Ultraviolet Luminosity Function: Stepwise Maximum Likelihood Method 


Ml 500 

V> (- ~ 4) 

P (- ~ 5) 
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H> (z « 7) 
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Note. — Magnitude bins with zero objects are shown as upper limits, calculated as 1/V e ///AM in that magni- 
tude bin for that redshift. 


density in each magnitude bin independently. 

For each magnitude bin and for each field, the likelihood 
was calculated (using Equations 5 and 6 above) that a given 
randomly drawn value of p(M) will give the observed num- 
ber of galaxies. The actual recorded value of <p(M) is that 
which maximizes the likelihood. While in practice, this yields 
very similar results as one would get by simply taking the 
observed number and dividing by the effective volume (con- 
sistent within a few percent for bins with more than a few 
galaxies), our approach has two advantages. First, in the 
limits where numbers are small, this approach is more accu- 
rate in that it properly accounts for the Poissonian likelihood. 
Secondly, this approach generates a full probability distribu- 
tion for the number densities in each magnitude bin, allow- 
ing for the derivation of accurate asymmetric uncertainties. 
Our SWML luminosity function determinations and best-fit 
Schechter functions are given in Table 5 and shown in Fig- 
ure 10. 

6. LUMINOSITY FUNCTION INTERPRETATION 
6.1. Evolution 

As shown in Figure 10, the qualitative shape of the SWML 
luminosity functions at all redshifts we consider here are sim- 
ilar, in that bright galaxies are rare and faint galaxies are rela- 
tively common. Additionally, when examining the Schechter 
fits (solid line), we see that they are consistent with the 
SWML determinations. The best-fit Schechter function pa- 
rameters (Table 4) surprisingly show little evolution in My V . 
However, from z = 4 to 8, the uncertainty on gets pro- 
gressively larger, to 0.4 (0.9) mag at z = 7 (8). This is easy 
to understand, as at all redshifts, our dataset contains galax- 
ies in only 1-2 bins brightward of M(j v . Ideally, one would 
prefer to have multiple bins in magnitude on either side of 
/V7j* jv to obtain robust constraints. As shown here, that will 
require a larger volume than we consider in this analysis. In 
Figure 1 1, we also fit the evolution of /V/ ( " v with redshift with 
a linear function, using our results at z = 4-8. We find that 
dM* /dz = -0.12 ±0.09; thus, our data do not support a sig- 
nificant evolution of M^ v with redshift. 

We also fit similar functions to see if we detect evolution 


in a and p*. As shown in Figure 11, we do see significant 
evolution in the faint-end slope a, with it becoming steeper at 
higher redshift, as da/dz = - 0.19 ±0.04 (4.8er significance). 
We see a similar significance in the evolution of the character- 
istic number density ip*, which evolves as dlogtp* / dz = - 0.31 
±0.07 (4.4cr significance). Thus, while My V does not signifi- 
cantly evolve with redshift from z = 4 to 8, both a and p* do, 
in that the number density decreases and the faint-end slope 
becomes steeper with increasing redshift. In particular, this 
decline in characteristic number density is by a factor of ~20, 
over a period of time of less than 1 Gyr. Although the steep- 
ening of the faint-end is consistent with previous studies (e.g., 
Bouwens et al. 2012), the un-evolving Mfjy and strong num- 
ber density evolution are the opposite of the picture presented 
in the literature just one year ago (e.g., Bouwens et al. 2007, 
2011b; McLure et al. 2013). This updated evolutionary pic- 
ture will be crucial when projecting number counts for future 
II ST and James Webb Space Telescope surveys. 

In Figure 12, we show both determinations of the luminos- 
ity functions together at all five redshifts, along with the joint 
confidence contours on v , a and p* . It is apparent that 
there is significant evolution in the luminosity function, with 
a drop in number density from z = 4 to 8, as well as a gradual 
steepening of the faint-end slope. The position of the “knee” 
of the luminosity function does not appear to evolve much, 
consistent with the results above that much of the evolution is 
in number density and not in magnitude. 

6.2. Impact of Magnitude Uncertainties 

By definition, our method of computing the luminosity 
function is dependent on magnitude binning, as we compare 
the observed number to that expected based on a given model 
in magnitude bins of width 0.5 mag. While galaxies close 
to one side of a magnitude bin have the potential to scatter to 
another bin, the typical uncertainties on the UV absolute mag- 
nitudes of galaxies in our sample are ^0.2 mag. Additionally, 
galaxies can shift both ways, thus while one galaxy moves out 
of a bin, another may move in, though this effect will not be 
symmetric given the shape of the luminosity function. In our 
results above, we had assumed that magnitude scatter does not 
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Fig. 11. — The evolution of the Schechter function parameters. We fit for evolution with redshift using our results from z = 4 to 8, and showing those from 
Bouwens et al. (2014) for comparison. Contrary to previous studies, we find no significant evolution in v . We find significant (>4cr) evolution in a, from 
steeper slopes at higher redshift to shallower slopes at lower redshift, and in the characteristic number density ip * , which evolves to higher values by a factor of 
20 x from z = 8 to 4. 



Fig. 12. — Left) The evolution of the UV luminosity function from z = 4 to 8, where the circles and lines denote our step-wise and Schechter-parameterized 
luminosity functions, respectively. Center) Contours of covariance between a and MJ V at z = 4, 5, 6, 7 and 8. The contours denote the 68% and 95% confidence 
levels, while the small circles show the best-fitting values. Right) Contours of covariance between p* and M at z = 4, 5, 6, 7 and 8. 


significantly impact our results. 

To investigate the impact of this assumption, we preformed 
another iteration of MCMC fitting to our data, here allowing 
galaxies to scatter between magnitude bins. At each step in 
the MCMC chain, a new value of M \ son was drawn for each 
galaxy from the 100 SED-fitting Monte Carlo simulation re- 
sults. The spread in these values encompassed both the pho- 
tometric scatter in the observed filters and the uncertainty in 
the photometric redshift (see §3.7). To compare to our fidu- 
cial luminosity function values, we recorded both the median 
Schechter parameter results and the median number density 
in each magnitude bin, as, unlike our fiducial MCMC run, 
these varied during each step as the magnitudes changed. At 
all redshifts, our fiducial values of the step-wise luminosity 
functions are consistent with these “magnitude-scatter” val- 
ues within 10% at M > -21.5, and typically within 2-3%. The 
sole exception is in the brightest bin (-22 atz = 4-6, and -21.5 
at z = 7 and 8), where our fiducial number density values are 
higher by ~ 15-20% (60% at z = 7, where there is only a sin- 
gle galaxy in this bin). We examined the Schechter fit, to see 
whether this bright-end difference affects our results. Values 
of both Myy and a derived when allowing galaxies to shift be- 
tween bins are consistent with our fiducial values within 0.1 
mag and <3%, respectively. We conclude that the relatively 
small ( ~20%) uncertainties in the absolute magnitudes of our 


galaxies do not have a significant impact on our luminosity 
function results. 

6.3. Non-parametric Evolution 

Given that our results show that the Schechter functional 
parameters may not be a robust method of tracking galaxy 
evolution (e.g., a non-evolving value of Mj*} v does not mean 
that the galaxy populations are not evolving), we examine the 
evolution in a non-parametric way. In Figure 13 we show the 
evolution of the step-wise luminosity function, plotting the 
number density corresponding to galaxies at Myv = -21 and 
-19 versus redshift. Fromz= 8 to 4 the abundance ofbrighter 
galaxies increases faster than faint-galaxies. This trend halts 
at z = 4, where bright galaxies have an approximately con- 
stant abundance down to z — 2, and then turns over. Faint 
galaxies, however, continue increasing in abundance down to 
z = 2, where they also turn over. This figure highlights the 
phenomenon of downsizing, where bright/large galaxies grow 
faster at early times (e.g., Cowie et al. 1996, see also Lund- 
gren et al. 2014). This is different from the expectation one 
would get simply from examining Schechter fits, as the lumi- 
nosity functions don’t evolve much over the range 2 < z < 
4 (e.g., Reddy & Steidel 2009). Given that the trends here 
mimic the evolution of the cosmic SFR density, we fit the 
function provided by Madau & Dickinson (2014) to our data 
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Fig. 13. — The number densities of bright (Myy = -21) and faint (Muy = 
-19) galaxies at a variety of redshifts. Our data are shown as large circles, 
and we also show results at high-redshift from Bouwens et al. (2014), Oesch 
et al. (2013) and Reddy & Steidel (2009). At low redshift, we show results 
as small circles from Amouts et al. (2005), Oesch et al. (2010b) and Cucciati 
et al. (2012). We fit the trend of c p with redshift, separately for our two 
magnitude bins, with the function given in Equation 7. The shaded regions 
show the 68% confidence ranges for each of the fits. The value of the slope of 
this function at high-redshift is significantly steeper for bright galaxies than 
for faint galaxies, showing that from z = 8 to 4, bright galaxies become more 
abundant at a faster rate than faint galaxies. This trend reverses at z = 4, 
where bright galaxies stop becoming more common. Bright galaxies peak in 
number density at z = 3. 1-3.2, sooner than faint galaxies, which peak at z = 
2.4-2. 7 (68% C. L.). At z < 2, the abundances of both populations plummet, 
in line with the evolution of the cosmic star- formation rate density. 


for both number densities, given by 


(1 +z) a 

VV)= A 1 . r/ I Twov-, ma S _1 M P C ~ 3 - 


The evolution with redshift is thus proportional to (1 +z)“ at 
low redshift, and (1 +z)“~ 7 at high redshift. Fitting the data in 
this way, we confirm that at z > 3, bright galaxies change in 
abundance faster, as (1 +z)~ 4 9±0,4 , than faint galaxies, which 
go as (1 +z) _3 - 3±0 - 3 . 

Another interesting aspect is to compare the trends ob- 
served to our predicted abundance of bright z = 9 galaxies 
(see §9). The trend observed here slightly overestimates our 
predicted z= 9 abundance, though if we assume the uncertain- 
ties on our z = 8 number density applies to z = 9, our trend is 
consistent with this prediction, in any case, this trend of abun- 
dance with redshift lends more weight to our expectation of a 
significant abundance ofbrightz= 9 galaxies. This figure ne- 
glects the impact of dust attenuation, as we are only looking at 
the observed UV magnitudes. The dust attenuation appears to 
be luminosity dependent (bright galaxies are dustier than faint 
galaxies, e.g., Bouwens et al. 2013), as well as being higher 
at lower redshift. Thus, correcting for dust would not only in- 
crease the abundance of bright galaxies more than that of faint 
galaxies, it would increase it by more at lower redshift, thus 
enhancing the differences between faint and bright galaxies at 
z > 4. 


6.4. Comparison to Previous Results 

Our result of a similarly bright value of M^ v at z = 6, 7 
and 8 is a dramatic change from previously published results. 
In Figure 10, we show the step-wise luminosity function re- 
sults from several relevant studies from the literature. Fig- 
ure 14 shows our uncertainty results, highlighting both the 
distribution of ) v and a from the MCMC chains, as well as 


the covariance between the two parameters, along with pre- 
vious determinations of M [ ) v and a (Bouwens et al. 2007; 
McLure et al. 2009; Ouchi et al. 2009; van der Burg et al. 
2010; Bouwens et al. 2011b; McLure et al. 2013; Schenker 
et al. 2013; Willott et al. 2013; Bowler et al. 2014). In this sub- 
section we compare solely to previous work - we reserve the 
comparison to the contemporaneous work by Bouwens et al. 
(2014) to §6.4.1 below. 

At z = 4 and z = 5, both our binned luminosity function 
data points as well as our Schechter function parameters are 
in excellent agreement with the ground-based study of van der 
Burg et al. (2010). We are also in excellent agreement with 
the ground-based study of McLure et al. (2009) at z = 5. We 
find good agreement with the space-based study of Bouwens 
et al. (2007) at z = 5, but at z = 4 the Bouwens et al. (2007) 
result lies outside our 2er confidence region on the Schechter 
function parameters, in that we prefer a shallower faint-end 
slope and a fainter value for 

At z = 6, our binned luminosity function data points are 
consistent within 1 -2a with the Bouwens et al. (2007) results 
at the faint end. At the bright-end, our data are higher than 
those from both ground-based studies (though again, typically 
only different at the l-2er level). This is somewhat counter- 
intuitive, as one may expect the ground-based studies to suffer 
a higher contamination rate, particularly for relatively fainter 
sources at higher redshift, due to their inability to resolve stars 
from galaxies, but it may also be explained due to an agressive 
sample selection required to minimize contamination. In any 
case, the differences are not highly significant, with the ex- 
ception of the brightest data point from Willott et al. (2013), 
which gives a number density at M = -22.5 of 2.7 x 10~ 8 
Mpc ’. While this is consistent with our upper limit at that 
magnitude, it is a factor of ^250 lower than our number den- 
sity only 0.5 mag fainter at M = -22 (see Table 5). Given the 
results at similar magnitudes at lower redshifts, it is highly 
unlikely that there is such a steep drop in number density 
over only a 0.5 magnitude interval, though future large area 
studies can better investigate the difference (Bowler et al. in 
prep). The larger discrepancy comes when comparing the 
Schechter function parameters. Specifically, Bouwens et al. 
(2007) find M* = -20.29 ± 0.19, and McLure et al. (2009) 
find M* =-20.04 ± 0. 12. Both values are significantly (2-3er) 
fainter than our derivation of M* = — 2 1 . 1 3^o |f . For the space- 
based study of Bouwens et al. (2007), this is understandable, 
as at that time only optical data were available, thus z = 6 
galaxies were selected via detections in only one band, and 
a robust determination of their UV absolute magnitudes was 
difficult. For the ground-based study of McLure et al. (2009), 
a cause for the difference is less clear, though certainly the 
different data being used plays a role. 

Comparing to previous works at z = 7, we find broadly sim- 
ilar results, in that our results are consistent with the derived 
number densities from previous studies, yet our Schechter fit 
prefers a much brighter value of My V . This is easier to under- 
stand, as a number of previous studies had less data available, 
and thus, utilizing smaller volumes, were unable to constrain 
the bright end (e.g., Bouwens et al. 2011b; Schenker et al. 
2013). The exception is the recent work by McLure et al. 
(2013), which used a similar volume as our study, though 
they used the CANDELS UDS field in place of our use of the 
CANDELS GOODS-N field. Examining the brightest data 
point from McLure et al. (2013) at M = —21, the number den- 
sity is about a factor of two below our data point. Flowever, 
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Fig. 14. — Confidence contours on our measured value of the faint-end slope a and the characteristic magnitude My V at z = 4, 5, 6, 7 and 8, with the light 
and dark shaded regions denoting 68% and 95% confidence. The large red circles represent our fiducial best-fit luminosity function parameters, while the other 
colored symbols denote results from previous studies, using the same symbols as in Figure 10 (with the addition of the results from Grazian et al. (2012), shown 
as the yellow triangle in the z = 7 panel, who fit a keeping fixed to -20. 14). In the z = 8 panel, we also show our best-fit result when fixing a > -2.3 as the 
white-filled red circle. The histograms to the top and side of each contour plot show the number of MCMC steps when a given value of M^ v or a was recorded, 
with the median value shown by the blue line. 


the discrepancy is mitigated by two factors. First, as dis- 
cussed by Bouwens et al. (2014), the use of fixed diameter 
circular apertures by McLure et al. (2013) systematically un- 
derestimates the fluxes for bright, more extended, galaxies. 
Bouwens et al. (2014) estimate the amplitude of this effect to 
be ~0.25 mag. Shifting the brightest McLure et al. (2013) 
data point by 0.25 mag brings it into agreement with our re- 
sults. Secondly, the CANDELS GOODS-N field appears to 
have an overdensity ofz= 7 galaxies. Specifically, when com- 
paring the number density of z = 7 galaxies in the GOODS-N 
Deep field to the GOODS-S Deep field in Figure 16, GOODS- 
N has a higher number density at all magnitudes. While we 
have not selected galaxy samples in the UDS, we can exam- 
ine this further by recomputing our z = 7 stepwise luminosity 
function using only the GOODS-S and FlUDF fields. At mag- 
nitudes fainter than -21 the results do not change appreciably, 
as the GOODS-N Deep and Wide fields lie on either side of 
our Schechter fit at those magnitudes. However, the results 
using only GOODS-S provide a number density ~33% lower 
at M= -21.5 than our fiducial luminosity function. This dif- 
ference is at the la level due to the large Poisson noise con- 
tribution in this bin, and thus is not highly significant. 

We also compare to several ground-based studies at z = 
7. Ouchi et al. (2009) identified 22 bright z ~ 7 candidate 
galaxies over ~^0.4 deg 2 . Their data points based on detected 


galaxies are consistent with our own, though their strict upper 
limits at M ~ -22 push their Schechter fit to a fainter value 
of My V = -20.1, although the large uncertainty (0.76) leaves 
My v consistent with our fit at only slightly more than the lcr 
level. The stepwise luminosity function from Castellano et al. 
(2010) based on deep HAWK-I data agrees well with our re- 
sults, while the results from the zFourGE medium band sur- 
vey of Tilvi et al. (2013) agree at M=-21.5, but differ by ~2 a 
at M= -20.5. 

Recently, Bowler et al. (2014) have made a significant im- 
provement in search volume from the ground, discovering 
34 luminous z ~ 7 galaxy candidates over 1.65 deg 2 , from 
the Ultra VISTA survey data over the COSMOS field (Mc- 
Cracken et al. 2012) and the UK1DSS survey over the UDS 
field (Lawrence et al. 2007). Broadly speaking, they are 
consistent with our results, and they are highly inconsistent 
with the previous determinations of MJ} V ~ -20 (Figure 10). 
There is a mild tension at M = -21.75, where the value of our 
Schechter fit at that point is 2cr higher than their derived num- 
ber density. However, this is their faintest magnitude bin, and 
is only ~50% complete, thus this data point relies the most on 
the completeness correction. In any case, the fact that Bowler 
et al. (2014) found z ~ 7 candidates out to very bright mag- 
nitudes gives us confidence that our brighter determination of 
My y is not necessarily dominated by cosmic variance in our 
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fields, but is a true feature of the z = 7 universe. However, our 
present uncertainty on My V of —0.4 mag makes it apparent 
that more data is needed to constrain this parameter further. 

At z= 8, we again find consistent number densities with pre- 
vious studies, though our larger volume allows us to find more 
rare, bright (M=- 21.5) galaxies than observed in some pre- 
vious surveys, pushing them to lower values of My V (though 
again here our uncertainty on M," v is large, so the difference 
in our determination is not significantly different from pre- 
vious studies). As noted above, in our fit of the z = 8 lumi- 
nosity function, we constrained My V to be fainter than -23, 
to avoid un-physically bright values, which tended to be pre- 
ferred in an unconstrained fit. We note two important points 
when comparing to previous studies. First, while our study 
did not utilize the pure parallel BoRG (Trenti et al. 2011) and 
HIPPIES (Yan et al. 2011) programs, our bright end is con- 
sistent with that from Schmidt et al. (2014), based on a de- 
termination of the z = 8 luminosity function over 350 arcmin 2 
of pure parallel data (for comparison, our search area at z = 8 
comprised —300 arcmin 2 ; Table 1). The multiple sight-lines 
of BoRG and HIPPIES leave their results less susceptible to 
cosmic variance effects, so the agreement implies that cosmic 
variance may not be strongly affecting our bright end, though 
we explore this in §6.5. 

A potentially larger difference between our results and 
those of previous studies is also seen at the faint-end, in that 
our faint-end slope is possibly steeper than previously found. 
However, our uncertainty is large, such that our result of 
a = -2.36^q4q is consistent with previous results of a « -2 
(e.g., Bouwens et al. 2011b; McLure et al. 2013; Schenker 
et al. 2013). Previous studies use galaxies as faint as M 
= -17.5 in their determination of the faint-end slope at z = 8. 
As discussed above, and shown in Figure 8, we find that we 
fall below 50% completeness at M > -18.5, thus we do not 
use galaxies fainter than that in our luminosity function de- 
terminations. While robust estimates of the number densities 
of galaxies at —18.5 < M < -17.5 would certainly improve 
the confidence on the faint-end slope, we use the same deep 
datasets as the other referenced studies (HUDF). We would 
expect the incompleteness to be similar between all studies, 
though it does depend on sample selection and the exact de- 
tails of the incompleteness simulations. In any case, con- 
straints on the faint-end slope at z = 8 should improve in the 
near future with further data from the Hubble Frontier Fields 
program. Our inclusion of the Hubble Frontier Fields paral- 
lel imaging, even though contributing only a small number of 
galaxies at z > 7, did improve the fractional error on the faint 
end slope by 2.3% and 7.8% at z = 7 andz= 8, respectively. 

Finally, there have also been theoretical estimates of the lu- 
minosity functions at these redshifts, most prominently from 
Jaacks et al. (2012a), who made predictions in good agree- 
ment with our observed luminosity functions. Specifically, 
their simulations also predict bright values ofM^y, of-21.15, 
-20.82 and -21.00 at z= 6, 7 and 8, respectively. They also 
found quite steep faint-end slopes, of -2.15^^, -2.30^Q 2 g 
and -2.5 Hq 22 at z = 6, 7 and 8, respectively. Within the un- 
certainties, these faint-end slopes are consistent with our own, 
though the apparent agreement at z = 8 is tantalizing (though, 
as mentioned above, we cannot constrain the slope to be so 
steep). Steep faint-end slopes of a ~ - 2 at these redshifts 
were also seen by Salvaterra et al. (2011) and Dayal et al. 
(2013), though both studies also predict a brightening in My V 
towards lower redshift which is no longer observed. 



Fig. 1 5 . — Left) Color-selection for z > 6.5 galaxies used by Bouwens et al. 
(2014) in the COSMOS, EGS and UDS fields, where HST 7-band data is not 
available. Right) Improved color selection in these fields with the addition of 
hypothetical 7 -band data. Of particular note is that without 7-band data, the 
z > 6.5 selection is potentially dominated by M, L and T dwarf stars. These 
clear out of the selection box with the addition of 7-band data. Additionally, 
galaxies with z < 6 spectroscopic redshifts from the literature (yellow boxes) 
move farther from the selection box, and are less likely to scatter in, with the 
addition of the 7-band data. 


6.4. 1 . Comparison to Bouwens et al. 2014 

Recently Bouwens et al. (2014) submitted a similar study 
of the evolution of the UV luminosity function at 4 < z < 10. 
Their sample of galaxies is larger than ours, as in addition 
to the datasets we use, they selected galaxies from the CAN- 
DELS COSMOS, EGS and UDS fields (though they did not 
use the HFF parallel fields). The data in these other CAN- 
DELS fields have a depth similar to the GOODS-S and N 
Wide fields, and thus are most useful for constraining the 
bright-end of the luminosity function. Comparing our results, 
while the agreement at z= 5 is excellent, the Bouwens et al. 
(2014) data points at z = 4 lie at higher number densities than 
our own for all but the brightest bins. These differences re- 
sult in a slightly steeper value of a at z= 4 (asouwens = -1-64 
versus ajhisSiudy = -1.56), but a significantly brighter value of 
M* y (-21 .07 versus -20.73). 

Bouwens et al. (2014) also selected galaxy samples at z = 

6, 7 and 8. Broadly speaking, they found similar results as we 
do at z = 6 and 7, in that previous studies determined values of 
My v which were too faint (Figure 14). However, investigating 
the actual data points in Figure 10, one can see that atz = 6 and 

7, the Bouwens et al. (2014) data points frequently he above 
our own. This is most significant in the brightest bin of their 
z = 7 luminosity function, which is 1 .7er higher than our point 
(interpolating amongst our M 1500 =-22 and -21.5 bins to de- 
rive a number density at their brightest magnitude of-21.86; 
their data point is 2.2cr higher if we compare it directly to our 
ATi 5 oo = -22 data point). At z = 8, the Bouwens et al. (2014) 
data points are more consistent with our own. However, they 
find both a fainter value of My v and a shallower faint-end 
slope. This is primarily due to their faintest data point, which, 
at M = -17.5, is well below our 50% completeness limit, and 
lies below the extrapolation of our measured luminosity func- 
tion, pushing them to a shallower slope. However, these dif- 
ferences at z = 8 are not significant, as Figure 14 shows that 
the Bouwens et al. (2014) results lie comfortably within our 
68% confidence contour on a and My v . If we constrained a at 
z = 8 during our fitting to be > -2.3, we obtain best-fit results 
similar to Bouwens et al. (2014, Figure 14). However, given 
the data at hand, there is no robust justification for such a con- 
straint, thus we do not include this in our fiducial luminosity 
function fits. 

The three CANDELS Wide fields used only by Bouwens 
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et al. (2014) lack space-based 7-band data, with HST data 
present in only four filters (V^o6, hu, J\is and These 

fields have deep ground-based optical data, although with 
much poorer angular resolution, and occasionally shallower 
depth than available with HST. Of particular worry is con- 
tamination by stars and/or brown dwarfs in these samples. 
The left panel of Figure 15 shows the color-selection plane 
for galaxies at z > 6.5 used by Bouwens et al. (2014) in the 
CANDELS COSMOS, EGS and UDS fields. While the selec- 
tion space used does include the likely colors of true z > 6.5 
galaxies, it also contains the bulk of M, L and T-dwarf tem- 
plate colors. As shown in the right-hand panel, by adding a 
single HST filter, the WFC3 7ios-band, stellar contaminants 
move out of the selection box, and lower-redshift galaxies 
move even further from the selection box. To mitigate stel- 
lar contamination, Bouwens et al. (2014) used both ground- 
based 7-band data, and the Source Extractor stellarity mea- 
surement. However, the ground-based data are presently not 
very deep, with Bouwens et al. (20 14) typically only detecting 
sources with 7 < 25.5 ( Muv, ^ < -21.5; see §3.4). Addition- 
ally, the stellarity measurement can only robustly distinguish 
point-sources from galaxies much brighter than the detection 
limit. Our test with the CANDELS //|f,o-band imaging in the 
COSMOS and EGS fields show that a robustly identified stel- 
lar sequence in the stellarity measurement is only possible at 
i/i6o < 25 (M ( /j / z= 7 < -22). In light of the apparent over- 
abundance of bright galaxies in the Bouwens et al. (2014) 
z = 6 and 7 samples compared to our results, we conclude 
that the higher quality data in our fields yield more robust and 
contamination-free measurements of the number densities of 
bright galaxies in the distant universe. 

6.4.2. Previously Published Measurement Uncertainties 

The differences in results, particularly on the characteristic 
magnitude M^ Y between our current study and previous stud- 
ies in the literature, are surprising, as in some cases the differ- 
ences are larger than what would have been expected given 
previously published uncertainties. In particular, Bouwens 
et al. (2011b) initially derived M* = -20.14 ±0.26 at z = 7. 
However, the data from Bouwens et al. (2011b) (Figure 10; 
green triangles), seem insufficient to constrain the bright end 
to such a relatively high precision. In particular, the evi- 
dence for a Schechter-like exponential decline at the bright 
end does not appear to be present from these early data, which 
is not surprising as this study was based on data only from 
the HUDF09 and ERS surveys, which is <20% of the vol- 
ume considered in our current work. To investigate this fur- 
ther, we performed another luminosity function fit to our data, 
using only data from the HUDF09 and ERS fields, finding 
M* v = -20.64 ± 0.92 at z = 7. Thus, without the CAN- 
DELS data, we find a somewhat fainter value for al- 
though fully consistent with our fiducial value, as well as the 
earlier value from Bouwens et al. (2011b). Likewise, at z = 
8, we find M(jy =-19.76 ± 1.40; fainter, but consistent with 
our fiducial z = 8 estimate. The same is not true for the work 
of Bouwens et al. (2014), who now find M* =-21.04 ±0.26 
at z = 7, which is brighter than their previous determination 
by ~2.5<r. Understanding the differences in the uncertainty 
computations between these studies is beyond the scope of 
our work, but we note that our current MCMC implemen- 
tation was designed to produce optimal uncertainties on the 
Schechter function fit parameters. As shown in Figure 11 
our current Schechter fit uncertainties are larger than those of 
Bouwens et al. (2014). While some of these differences may 


be due to the fact that they used a larger volume (including all 
five CANDELS fields), the different methods of computing 
the uncertainties likely play a role. 

6.5. Cosmic Variance 

The impact of cosmic variance on our measurement of the 
luminosity function is minimized due to our use of several 
fields, which are split into four widely separated regions of 
the sky. However, as shown in Figure 16, there is significant 
variance between the different fields, particularly at z > 5. To 
estimate the effect of cosmic variance on our derived num- 
ber densities, we used the QUICKCV calculator provided by 
Newman & Davis (2002). For a given survey geometry, this 
program returns the fractional error in a count due to cosmic 
variance. For this calculation, we estimated the fractional er- 
ror separately for GOODS-S, GOODS-N, MACS-0416 par- 
allel, and Abell 2744 parallel fields, adding the variances in 
quadrature to derive a final value of ocv for a given redshift 
bin. In the GOODS-S field, we include the area from the three 
HUDF09 fields, as even the parallel fields are separated by 
only a few arcmin from the GOODS-S proper. For the input 
survey geometries, we estimate rectangular regions of the ap- 
proximate shape of the GOODS fields, with an enclosed area 
equal to the GOODS-S Deep+Wide+ERS+HUDF09 fields for 
GOODS-S, and the GOODS-N Deep+Wide for GOODS-N. 
The field geometries were thus 10.2' x 15.03' for GOODS-S, 
9.51' x 14.65' for GOODS-N, and 2.1' x 2.1' for each of the 
HFF parallel fields. With these inputs, we find values of ocv 
of 0.111, 0.106, 0.115, 0.124 and 0.133 atz= 4, 5, 6, 7 and 
8. The impact of cosmic variance on our sample is thus at the 
~10-13% level. 

To assess the impact of cosmic variance on our measured 
luminosity functions, we compare this uncertainty to the Pois- 
son noise from our step-wise luminosity functions. At all red- 
shifts, the data at M < -21 are dominated by Poisson noise 
(M < - 20.5 at z ± 6), thus we do not expect cosmic vari- 
ance to be dominating the uncertainties on the bright-end of 
the luminosity functions derived here. However, cosmic vari- 
ance may play some role at the faint-end, where we are re- 
stricted to small fields. At z = 7, there does appear to be a 
step in the stepwise luminosity function at M > -18.5, where 
the M = —19 point is below our best-fit Schechter function, 
and the M = -18.5 point is above. At M > -18.5 our data 
come from only the HUDF main field, thus this break rep- 
resents the point where we become reliant on a single small 
field. QUICKCV estimates that a single HUDF-sized field at 
z = 7 has a cosmic variance uncertainty of 36.2%. Comparing 
to our Poisson noise estimate of 30% uncertainty at M = — 18.5 
atz = 7, cosmic variance may bear some responsibility for this 
discontinuity in the luminosity function at the faint end. Fu- 
ture measures of the luminosity function at M > -18.5 from 
the Hubble Frontier Field lensing program may improve these 
constraints, but while faint galaxies may be found, the vol- 
umes will still be incredibly small. Thus, robust constraints 
on the number densities at this faint level at z > 7 may need 
to wait until the James Webb Space Telescope. 

We note that the cosmic variance estimates from 
QUICKCV are for dark matter only, and thus assume a bias of 
unity. The bright galaxies we observe are likely more biased, 
thus these estimates are lower limits of the impact of cosmic 
variance on our results. While multiple fields could in prin- 
ciple allow one to empirically measure the effect of cosmic 
variance, we have only two independent large fields, thus this 
is not robustly possible with our current dataset. 
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Fig. 16. — The rest-frame UV luminosity functions at each redshift for each 
Upper limits are shown for magnitude bins with zero detected galaxies. 

6.6. Do the Data Support a Schechter Function? 

When allowed to choose any value of M^ v , our z = 8 
Schechter function fit preferred very bright values of M,j V , 
such that all observed data points lay on the faint-end slope 
part of the function. This implies that the z = 8 luminosity 
function may be consistent with a single power-law. Such a 
functional form is what one might expect when the feedback 
effects which govern the bright end at lower redshift (mainly 
feedback due to accreting supermassive black holes) disap- 
pear, or if dust attenuations ceases to be a factor. Bowler et al. 
(2014) recently postulated that the z= 7 luminosity function 
is better fit by a double-power law, rather than a Schechter 
form. At z = 7, our step-wise data appears consistent with the 
Schechter fit out to the brightest magnitudes we cover. To see 
whether our data show a preference for a Schechter functional 
form at all redshifts, we performed three fits to the data - a 
Schechter fit, a single power law, and a double power law. To 
place these fits on equal ground, we found the best-fit param- 
eters for each function using a simple maximum likelihood 
routine. For the Schechter fit, we use the function shown in 
Equation 4, investigating a range of M^ v with AM= 0. 1 mag, 
and a with Aq = 0.02. We approximated a single power law 
using the Schechter functional form with fixed at -30. 
For the double power-law, we used the form given in Equation 
2 of Bowler et al. (2014), which is similar to the Schechter 
function at the faint-end, but replaces the bright end with a 
second power law with slope /?. In all cases, if* is found as 
the normalization such that the total number of expected ob- 
jects for a given function is equal to the number observed. The 
likelihood that a given functional form represents our data is 
calculated in an identical manner as in §5.1, using Equations 
5 and 6. 

To compare the results from these fits at each redshift, we 
used the Bayesian information criterion (BIC). This is simi- 
lar to a x 2 statistic, except that it takes into account both the 


subfield. The solid line denotes the best-fit Schechter function fit at each redshift. 

number of data points and the number of free parameters. For 
a model to be preferred over a competing model, it must have 
a BIC lower by at least 2. This is sensible, as adding a free pa- 
rameter must yield a better fit for that model to be preferred. 
The BIC is calculated as 

BIC = -2 \n(C) + k In (TV) (8) 

where N is the number of data points, and k is the number 
of free parameters (Liddle 2004). For the Schechter, double 
power law, and single power law fits, the number of free pa- 
rameters are 2 (M{j v , a), 3 (My V , a, (3) and 1 (a), respectively 
(we do not count ip* as a free parameter as it is a normaliza- 
tion). The number of data points is the number of galaxies in 
our sample used in the fit, which is restricted to those brighter 
than the 50% completeness limits discussed above. This gives 
N = 2788, 1812, 605,221 and 47 galaxies at z = 4, 5, 6, 7 and 
8, respectively. 

The results from this analysis are shown in Table 6. A dif- 
ference in the absolute value of the BIC of 2 is interpreted 
as “positive” evidence, while a difference of 6 or higher is 
“strong” evidence, both in favor of the model with the smaller 
value. In Table 6, in addition to the value of the BIC, we 
show the difference between the BIC values for the Schechter 
versus double power-law, and Schechter versus single power- 
law. In this formalism, a negative difference is in favor of the 
Schechter function. Comparing the Schechter function ver- 
sus the double power-law, we find that a Schechter form is 
strongly preferred to either a double or single power law at 
z = 4-7. This is not surprising, as there is clearly a deficit 
of observed galaxies at the bright end when compared to the 
best-fit power law (Figure 10). However, no such deficit is 
visible at z = 8, and this is confirmed as both the Schechter fit 
and the single-power law fit have effectively identical values 
of the BIC. We conclude that our data support an exponential 
decline at the bright end of the luminosity function at z = 4, 
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Fig. 17. — Rest-frame ultraviolet luminosity functions at a = 4—8, comparing our observations to the semi-analytic models (SAMs) of Somerville et al. (2012) 
(S12). The crosses mark the models of S12 both with and without dust. In order for the models to be consistent with the observations, it is clear that some dust is 
needed at all redshifts (except perhaps atz= 8), particularly at the bright end. Using a simple model for dust attenuation, our results suggest that the dust-to-metal 
ratio or dust geometry must change as a function of redshift, a continuation of a trend already pointed out in previous studies. 


TABLE 6 

Comparison of Luminosity Function Fits 


Redshift 

BIC 

Schechter 

BIC 

Double 

BIC 

Power 

ABIC 

Sch-Dou 

ABIC 

Sch-Pow 

4 

358 

377 

641 

-18.5 

-283 

5 

540 

562 

694 

-21.7 

-153 

6 

350 

361 

376 

-11.1 

-25.7 

7 

225 

234 

235 

-9.28 

-9.83 

8 

86.9 

92.3 

86.7 

-5.36 

0.26 


Note. — The comparison of the Bayesian information 
criterion statistic for fits to our z = 4, 5, 6, 7 and 8 luminos- 
ity functions using a Schechter, double power-law and single 
power law functional fonn. A difference in the absolute value 
of BIC between two models of >2 (6) is positive (strong) 
evidence for the preference of one model over another. A 
Schechter function is strongly preferred over a single power 
law at all redshifts except z = 8, where our data cannot distin- 
guish between the two models. 

5, 6 and 7. At z = 8, we do not see any evidence for a decline 
in the bright end, at least out to M = -21.5. Further data are 
needed to tighten these constraints, to show whether one can 
either detect, or rule out a decline at the bright end at z= 8. 
If the latter ends up the case, it could indicate either a signif- 
icant change in the halo masses of bright galaxies, a drop in 
dust attenuation in bright galaxies, or a change in the physics 
governing the feedback in bright galaxies in the distant uni- 
verse. 

7. COMPARISON TO SEMI-ANALYTIC MODEL PREDICTIONS 

In this section we compare our observations with predic- 
tions from theoretical models set within the predominant A 
Cold Dark Matter (ACDM) paradigm. All such models, 
whether based on numerical hydrodynamics or semi-analytic 


techniques, currently rely upon phenomenological "sub-grid” 
recipes to treat the physics on scales smaller than those that 
can be directly resolved. These processes include star for- 
mation and feedback from massive stars, supernovae, and su- 
permassive black holes. The phenomenological recipes are 
parameterized and must be empirically calibrated. Here, we 
compare our new observations at z = 4 - 8 with predictions 
from the models presented in Somerville et al. (2012, here- 
after S12). The sub-grid recipes in these models have been 
calibrated using a set of observations at z ~ 0, and Somerville 
et al. (2012) presented a comparison with available observa- 
tions fromz ~ 0-5. It is therefore very interesting to test these 
model predictions — with no re-tuning of the free parameters 
controlling physical processes — in the higher redshift regime 
probed in this work. 

Figure 17 shows our estimates of the rest-UV luminosity 
function compared with the S12 SAM predictions with and 
without dust. It is already interesting that the dust-free model 
predictions are even in plausible agreement with the observa- 
tions; i.e., the model predictions lie above the observations at 
all luminosities and redshifts. Next we can ask the question: 
what characteristics must the dust extinction have in order to 
be consistent with the observations? One can immediately 
see that the dust extinction must be differential with both lu- 
minosity (more luminous galaxies are more extinguished) and 
redshift (galaxies are less dusty at higher redshift). We use a 
simple approach to model the dust extinction: as in S12, we 
assume that the face-on dust optical depth in the F-band is 
given by r 0 y = r dust m co i d Z cold /r gas 2 , where m cold is the mass of 
cold gas in the disk, Z co i d is the metallicity of the cold gas, r gas 
is the exponential scale radius of the gaseous disk, and r dust is 
a normalization parameter. The values of m co id, Z co ] d , and r gas 
are predicted by the SAM. We treat r dust as a free parameter. 
We then assign random inclinations to our galaxies and use a 
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“slab” model to compute the inclination-dependent extinction 
(see S12 for details). We use a Calzetti attenuation curve to 
compute the attenuation at 1500 A. 

In SI 2, we showed that if we normalize r ( i llst to match ob- 
servations at z ~ 0 and use a fixed value, our model over- 
predicts the dust extinction at higher redshift. Similarly here, 
we find that the empirical redshift-dependent function for Td us t 
adopted in S12 based on observations at z < 5 overpredicts 
the extinction at z > 5. We empirically adjust Td us t to obtain 
a good fit to the observed LF in the five redshift bins shown, 
and find that Td us t oc exp(-z/ 2 ) produces a reasonably good 
fit over this redshift range, where z is redshift. This may be 
physically interpreted as either a changing dust-to-metal ra- 
tio, or a systematic evolution in the dust geometry relative to 
our simple slab model. The required luminosity and redshift 
dependence of the dust extinction is in qualitative agreement 
with observational conclusions drawn based on the UV colors 
(Finkelstein et al. 2012b; Bouwens et al. 2013). 

In future work, we plan to investigate whether the dust ex- 
tinction parameters derived from SED fitting on the observa- 
tions are consistent with the empirical SAM requirements. In 
addition, we plan to use these models, which plausibly match 
the observed UV luminosity functions, to make predictions 
for the clustering, stellar fractions, and other properties of 
high redshift galaxies. We will also show the results of vary- 
ing the sub-grid recipes for star formation and feedback, to 
illustrate what physical insights can be gained from these ob- 
servations. For the moment, however, it is intriguing that the 
models that were developed to explain galaxies at a very dif- 
ferent epoch are plausibly consistent with these new observa- 
tions. 

8. EVOLUTION OF THE COSMIC STAR-FORMATION RATE 
DENSITY 

While the evolution of the shape of the luminosity function 
can provide interesting constraints on the physics of galaxy 
evolution, the integral of the luminosity function provides a 
key measure of the total number of UV photons produced at a 
given redshift. This is a key constraint in two ways. First, the 
integral of the total SFR density provides a key check against 
the measured stellar mass density (e.g., Bouwens et al. 2013). 
Secondly, assuming a conversion between UV and ionizing 
photons, this measure can determine whether galaxies are pro- 
ducing enough ionizing photons to reionize the universe at 
a given redshift (e.g., Madau et al. 1996; Finkelstein et al. 
2012a; Robertson et al. 2013). 

We calculated the luminosity density at a given redshift, 
integrating down to Myy =-17. This is approximately the 
magnitude of the faintest galaxy in ourz = 8 sample, and also 
facilitates comparison with recent works which use a similar 
magnitude limit. Galaxies likely exist beyond this magnitude 
limit (e.g., Trenti et al. 2012; Alavi et al. 2014), which we will 
consider in the next subsection. We utilized the results of our 
MCMC luminosity function fitting chain to derive a robust es- 
timate of both the rest- frame U V specific luminosity density 
(puv, in units of erg s 1 Hz 1 Mpc ! ) and its uncertainties. In 
each step of the chain, we calculated puv by taking the lumi- 
nosity function from the best-fit Schechter function parame- 
ters for that step, and integrating it from -23 < M 1500 < -17. 
To convert this number to a SFR density, we use the rela- 
tion adapted from Kennicutt (1998, psfr = 1.25 x 10 - 28 puv), 
which converts the specific UV luminosity density to a SFR 
density ( psfr ), assuming a Salpeter IMF and a constant star- 
formation history over >100 Myr. The original coefficient 


TABLE 7 

Rest-Frame UV Luminosity Densities and SFR 
Densities 


Redshift 

log puv 

log SFR Density 


(ergs s -1 Hz -1 Mpc -3 ) 

(M 0 yr 

1 Mpc 3 ) 


Observed 

Observed 

Dust-corrected 

4 

26 26 +001 
zo,zo - 0 . 01 

1 cQ+0.01 
-i.~>?_o.oi 

-1 03 + °- 23 
1 -0.21 

5 

26 17+ 001 
zo - 1 -0.01 

-1 69 +001 
1,0 -0.01 

-i-20T8:i§ 

6 

9c 88+0-02 

ZD.OO-0 02 

1 Q7+0.02 
~ L - y '-0.02 


7 

9 c 77+O.O6 
ZJ - / - 0.06 

9 OQ+0 06 

-z.u ?_ 0 .06 


8 

25 65 +019 
-0.19 

-2 20 +019 
• -0.19 

-2-2(tS;g 


Note. — All values have been computed down to M uv = 

-17. The dust correction was derived based on the values 
of E(B-V) derived from SED fitting, with the dust-corrected 
SFR densities including an uncertainty term from the spread 
of extinction values at a given absolute magnitude. The SFRs 
were computed assuming the Kennicutt (1998) conversion 
from the UV luminosity density (puv), assuming a Salpeter 
IMF, and a constant star-forming population with age >100 
Myr. 

from Kennicutt (1998) was 1 .4; however, updated stellar pop- 
ulation models (e.g., Bruzual & Chariot 2003) imply a some- 
what smaller value. We chose a value of 1 .25 to be consistent 
with the assumptions used in Bouwens et al. (2014), though 
we note an even lower coefficient of 1.15 was used in Madau 
& Dickinson (2014). The quoted value of puv or psfr is the 
median of the values recorded from all of the MCMC steps, 
while the 68 % confidence range is taken to be the central 68 % 
of values. 

Although the UV luminosity is a relatively easy observable 
in this epoch, the major drawback in its use as a SFR indi- 
cator is its susceptibility to attenuation by dust. As a bevy 
of recent work has shown, this dust correction is important 
even out to z ~ 7-8 (e.g., Finkelstein et al. 2012b; Dunlop 
et al. 2013; Bouwens et al. 2013, see also §7). To calculate 
the total SFR density, we corrected the observed SFR density 
for extinction using a new iteration of SED fitting, including 
the deep Spitzer/IRAC data (§3.6), which is a crucial probe 
of the rest-frame optical light, providing better constraints on 
the dust attenuation. Using these updated extinction results, 
we calculated a sigma-clipped median and standard deviation 
for the best-fit extinction values at a given redshift in four 
magnitude bins: < - 21, -21 to -20, -20 to -19, and -19 
to -17. We recover previously observed trends that dust ex- 
tinction lessens with both increasing redshift and decreasing 
UV luminosity (e.g., Finkelstein et al. 2012b; Bouwens et al. 
2013). The values of E(B- V) for bright galaxies decreases 
from 0. 15 at z = 4 to 0.02 at z = 7, and for faint galaxies from 
0.06 atz= 4, to 0.0 at z= 7. The small numbers, limited wave- 
length coverage, and faint magnitudes of z = 8 galaxies make 
it difficult to measure their extinction, therefore we assumed 
E(B-V) = 0 for allz = 8 galaxies. The spread in E(B-V) val- 
ues at all redshifts and luminosities is ~ 0 . 1 , thus we assume 
this value in all cases (with the exception of z = 8 , where we 
fix E(B- V) to zero). To include this uncertainty in E(B— V) 
in our derived dust-corrected SFR density, in each step of the 
chain we draw a new value of E{B — V) for a given redshift 
and magnitude bin, modifying the fiducial value by a number 
drawn from a Gaussian distribution with a standard deviation 
equal to the E(B-V) spread of 0.1. The values of puv and 
the observed and dust-corrected values of psfr are given in 
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Fig. 18. — The evolution of the cosmic star-formation rate density, derived by integrating the best-fit Schechter function at all redshifts to Muv < - 17. Our 
data are shown as large circles. To extend this analysis to lower redshifts, we also show the values at z ~ 2-3 from Reddy & Steidel (2009), and from z = 0-2 
from Amouts et al. (2005). For both studies, we integrated the published best-fit Schechter function parameters to -17 to derive the uncorrected values of puv- 
We used the published ratio of the dust-corrected - to - unobscured values of puv from Reddy & Steidel (2009) to calculated the dust-corrected values at z ~ 
2-3. Atz < 2, we used the dust correction from Schiminovich et al. (2005), which assumes a constant value of Auy = 1.8 at all redshifts. We used Equation 7 
to fit the observed trends, deriving the uncertainties on this fit via 10 4 Monte Carlo simulations, showing as the shaded region the 68% confidence range from 
these fits. The total (dust-corrected) SFR density evolves with (l+z) -4 ' 3 ^ 0,5 fromz = 3-8. The green symbols show the high-redshift results from Bouwens et al. 
(2014) and Oesch et al. (2013), which were not included in the fit, but are consistent with the observed trends at the ~ I a level. 


Table 7. These values do not include potential sub-millimeter 
galaxies which lie below our rest-frame UV detection limits. 
However, as we are observing at z > 4, we expect their im- 
pact on the total cosmic SFR density to be minimal (see Table 
7 from Bouwens et al. (2012)). 

In Figure 18, we show our derived values of the cosmic 
SFR density. Our results are for the most part consistent with 
those of Bouwens et al. (2014), although our observed val- 
ues are lower by about a few a at z = 4 and 5, likely due to 
our shallower faint-end slopes at these redshifts. To study the 
evolution of puv with redshift, we fit the function provided by 
Madau & Dickinson (2014), given in Equation 7. Although 
we have included lower-redshift data in our fit, we do not dis- 
cuss here our results for the low-redshift slope or the peak red- 
shift, as these are better obtained from Madau & Dickinson 
(2014), who use a compilation of several sources, including 
far-infrared observations. As we are adding data at high red- 
shift, it is interesting to examine the trends there. We find that 
atz> 3 the uncorrected values ofpuv evolve as (l+z) -2 ' 4±0 - 3 , 
while the dust corrected values evolve as (1 -Fz) -4 3±0 5 . Most 
interesting, the observed trend of the evolution of the total 
SFR density is consistent within 1 a of the published results 
at z = 9 from Oesch et al. (2013) and at z= 10 from Bouwens 
et al. (2014). We thus do not see any evidence that there is a 
break in the evolution of the cosmic SFR density beyond z = 
8 . 

8.1. Constraints on Reionization 

Although it is presently generally assumed that galaxies 
dominated the ionizing photon budget for the reionization of 


the IGM, it has been difficult for observations to obtain robust 
proof. Analyses of the IGM via line-of-sight quasar observa- 
tions have been able to show that reionization was likely com- 
plete by z ~ 6 (e.g., Fan et al. (2006), though see Mesinger 
(2010) and Becker et al. (2014)). Additionally, observations 
of the cosmic microwave background (CMB) radiation con- 
strain the total optical depth due to electron scattering, which, 
while it cannot directly inform the duration of reionization, 
it can give an estimate of the reionization redshift ( z,. CTO „ ) if 
reionization was instantaneous. The results from the Wilkin- 
son Microwave Anisotropy Probe (WMAP) 9-year dataset 
give r es = 0.088 ± 0.014, which corresponds to z reion = 10.6 ± 
1.2 (Hinshaw et al. 2013). 

The primary reason for the current uncertainties in the con- 
tribution of galaxies to reionization lies in the uncertainty in 
the faint-end slope measurements, and also in the assump- 
tions of the escape fraction of ionizing photons (f esc ) and the 
clumping factor in the IGM (C). The clumping factor is pri- 
marily constrained theoretically, but most studies agree that 
it is low (<6) at high redshift (e.g., Faucher-Giguere et al. 
2008; Pawlik et al. 2009; McQuinn et al. 2011; Finlator et al. 
2012). To infer from observations of galaxies a number of es- 
caping ionizing photons, one first needs to take the observed 
UV light, and assume an IMF and a metallicity. Then, to 
calculate the number of these ionizing photons available for 
reionization, one then needs to multiply by an assumed value 
of f esc . It is difficult to constrain f esc directly at high-redshift, 
as the correction for intervening IGM absorption systems is 
extremely high at z > 4. Significant effort is being expended 
on observationally constraining / esc atz < 4. Although bright 


28 


Finkelstein et al. 



0.12 
0.10 
0.08 
0.06 ^ 
0.04 

0.02 
0.00 

4 6 8 10 12 14 

Redshift 



Fig. 19. — Left) The specific luminosity density (puv) versus redshift (similar to Figure 3 from Finkelstein et al. (2012a)). Here we show our luminosity 
functions integrated down to M< —13 as blue circles. The cyan circles denote the value of puv when we integrate down to our 50% completeness limit (—18 at 
z = 7). Recent results from Bouwens et al. (2014) at z ~ 10 are shown in green, with the upper and lower squares representing limiting magnitudes of —17 and 
-13, respectively. The wide gray curves denote the value of puv needed to sustain a fully reionized IGM at a given redshift, for a given ratio of the clumping 
factor C over the escape fraction of ionizing photons fesc (Madau et al. 1999). The thin blue curve shows our fiducial value of C =3 and f esc = 13%. Right) 
The volume ionized fraction, xhii, of the IGM which can be sustained given the integral of our luminosity functions at z = 4-8 (as well as that at z = 10.4 from 
Bouwens et al. (2014)). We assume the luminosity function extends to Muv = — 13, C=3 and f esc — 13% (this escape fraction is the highest that does not violate 
constraints set by the Lya forest at z = 6; Finkelstein et al. 2012b). We plot constraints on xhii from spectroscopy of quasars at z < 6 from Fan et al. (2006) and 
at z = 7 from Bolton et al. (201 1). The blue circle denotes constraints on xhii from the evolution in the Lya luminosity function from z = 5.7 to 6.6 from Ouchi 
et al. (2010), while the blue bar denotes the range of xhii values inferred from z = 7 follow-up Lya spectroscopic studies (e.g., Pentericci et al. 2014; Tilvi et al. 
2014). The instantaneous redshift for reionization from WMAP (10.6 =b 1.2) is indicated by the orange rectangle. The derived 50% and 90% xhii redshifts from 
the study of Kuhlen & Faucher-Giguere (2012) are shown in green. The right-hand axis corresponds to the hatched regions, which show the Thomson optical 
depth to electron scattering ( r es ) as predicted by our integrated luminosity functions (blue) compared to WMAP (orange). Compared to previous results, the 
improved constraints on the luminosity functions yield a tighter range of possible reionization histories. Broadly speaking, we find a picture where the universe 
is fully ionized by z = 6, with the neutral fraction becoming non-negligible at z > 7, with r es consistent within 1.3cr of the WMAP9 value. 


galaxies at z ~ 1 have very low escape fractions (relative to 
the UV emission) of f esc ,rei < 2% (Siana et al. 2010), escap- 
ing ionizing emission has been observed from small fractions 
of galaxies probed by studies at z ~ 3-4 (e.g., Steidel et al. 
2001; Shapley et al. 2006; Iwata et al. 2009; Vanzella et al. 
2010; Nestor et al. 2011), though some ground-based studies 
may suffer from contamination by intervening sources (e.g., 
Vanzella et al. 2012). Recent results imply that escape frac- 
tions from star- forming galaxies at z ~2— 3 range from 5-20%, 
with lower-mass galaxies, especially those with Lya in emis- 
sion, having a greater likelihood of having detectable escap- 
ing ionizing emission (e.g., Nestor et al. 2013; Mostardi et al. 
2013). 

Finkelstein et al. (2012a) used measurements of the emis- 
sion rate of ionizing photons from observations of the Lya 
forest in quasar spectra to place an upper limit on f esc from 
galaxies. Assuming that the rest-frame UV luminosity func- 
tion extended down to M\jy = -13, the escape fraction must 
be f esc < 13% to avoid violating the Lya forest measurements 
of Bolton & Haehnelt (2007). Using this value, and assum- 
ing C = 3, the luminosity functions available at the time were 
consistent with a wide range of reionization histories, includ- 
ing an end redshift as late as z < 5, and an ionized fraction 
at z ~ 7 from 30-100%. Kuhlen & Faucher-Giguere (2012) 
and Robertson et al. (2013) did similar analyses, folding in 
additional observables (e.g., the Lya forest and CMB), found 
that in order to complete reionization by z ~ 6 , the luminosity 
function must extend much deeper than can presently be ob- 
served, and/or the average escape fraction must be higher at 
higher redshift. 


Here, we use our updated luminosity functions to reexam- 
ine the contribution of galaxies to reionization. Figure 19 
shows both the observable specific UV luminosity density 
(Puv), which we define to be that above our 50% complete- 
ness limit, as well as the total puv, which we define as the 
integrated luminosity function down to M \ 500 = -13. We then 
compare these values to the critical number of UV photons 
necessary to sustain an ionized IGM at a given redshift, taken 
from Madau et al. (1999). This figure is similar to Figure 3 
from Finkelstein et al. (2012a), thus we refer the reader there 
for more details. Effectively, these critical curves depend on 
assumptions about the stellar IMF, metallicity, f esc and clump- 
ing factor. The first two are responsible for the conversion 
from observed UV photons to intrinsic ionizing photons. We 
assumed a Salpeter IMF, and the width of the curves denote 
the impact of changing the metallicity from 0.2 < Z/Zq < 1.0. 
We show several curves for the reader’s choice of the ratio of 
Cj f esc . Here, we use a fiducial value of C = 3 and f esc = 13%, 
consistent with Finkelstein et al. (2012a). 

The right panel of Figure 19 shows the ionization history 
of the IGM, comparing our derived value for the total specific 
UV luminosity density to our fiducial model ofC = 3 and f esc = 
13%, folding in the values at z= 10.4 from Bouwens et al. 
(2014) to extend our analysis beyond z = 8 . Our luminosity 
functions are consistent with a reionization history that starts 
at z ~ 11, and ends by z > 5. Although the exact value of 
the volume ionized fraction in the IGM is uncertain between 
these redshifts, due to the persistent uncertainty in the faint- 
end slope, our results imply the following constraints (given 
the caveat of our assumptions). At z = 6 , we can constrain 
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X//// > 0.85 (ler), while out to the limit of our observations 
at z = 8 the data are still consistent with a fully ionized IGM 
( 68 % C.L. of 0.15 < xhii =< 1.0). We find a midpoint of 
reionization (x Hn = 0.5) of 6.7 < z < 9.4 ( 68 % C.L.). 

Broadly speaking, measurements from quasar spectra as 
well as from Lya emission from galaxies support a reioniza- 
tion scenario consistent with what we derive (Figure 19). The 
constraints from Lyo- emission are heavily model dependent, 
and studies claiming a very low value of xhu may be assum- 
ing a velocity offset of Lya from systemic which is too high 
(e.g., Stark et al. 2014). The one measurement which is in 
slight tension with our results is that from WMAP9. From 
our fiducial reionization history, we find r es = 0.063 ± 0.013. 
As can be seen from the juxtaposition of the 68 % confidence 
regions in 19, the tension between the value of T es inferred 
from our results and that measured by WMAP9 is only slight, 
at 1.3cr. 

Future observations are necessary to improve the con- 
straints on reionization from galaxies. Specifically, more ro- 
bust measurements of the faint-end slope a at z = 6-8 can 
dramatically shrink the uncertainties on puv- subsequently re- 
ducing the width of our plausible values of xhii- Likewise, im- 
proving the measurements at z > 9 will inform us on whether 
the ionization fraction of the IGM at that early time was sig- 
nificantly non-zero. Even a small contribution (~10%) to xhii 
at early times will erase any discrepancy between our cur- 
rent observations and those from WMAP. The Hubble Fron- 
tier Fields program will improve both of these areas, though 
definitive results will likely not be obtained until the James 
Webb Space Telescope (JWST) era. 

9. EVOLUTION OF GALAXIES AT Z > 9 

Studies of galaxies at z > 9 are now only in their nascent 
phase, but HST surveys such as CANDELS and UDF12 are 
beginning to probe this early epoch. The first robust re- 
sults on galaxies in this epoch were published by Ellis et al. 
(2013), who used the new F140W data in the HUDF from 
the UDF12 program to discover galaxies at z ~ 9. This fil- 
ter allows z ~ 9 galaxies to be detected in two bands (F 1 40 W 
and F160W), dramatically reducing the contamination due to 
noise from F160W-only studies alone (§3.8.1; c.f, Bouwens 
et al. 2011a). Ellis et al. (2013) discovered the first robust 
sample at z > 8.5, finding seven candidate galaxies. McLure 
et al. (2013) followed this up with an analysis of the z = 9 
luminosity function, finding number densities at the faint-end 
(Mjjv ~ —18) only slightly lower than at z= 8 . Oesch et al. 
(2013) also analyzed the z = 9 luminosity function, also find- 
ing seven z ~ 9 candidate galaxies in the HUDF. Although 
the GOODS-S field lacks the F140W data necessary to de- 
tect potential z = 9 galaxies in two-bands, Oesch et al. (2013) 
added the full CANDELS/ERS GOODS-S field to improve 
their constraints at the bright end. However, they found no 
z = 9 candidates in this larger field. Their published luminos- 
ity function is consistent with that from McLure et al. (2013) 
at the faint end. Bolstered with their additional constraints 
due to the inclusion of the non-detections from the larger 
GOODS-S field, Oesch et al. (2013) fit a luminosity function 
(keeping the faint-end slope and normalization fixed), finding 
a surprisingly faint value for Mjjy of -18.8 ± 0.3. However, 
this derivation was based on the assumption that the luminos- 
ity function shows luminosity evolution at z > 6 — a trend 
which we have now shown to be unlikely. Given this new in- 
sight, as well as the presence of a plethora of bright galaxies 


at z = 7 and 8 , we consider it likely that the Oesch et al. (20 1 3) 
estimate of the bright-end of the z = 9 luminosity function is 
underestimated. 

A number of recent papers have described empirical evi- 
dence that galaxies at high redshift have star-formation histo- 
ries that increase with time (e.g., Papovich et al. 2011; Salmon 
et al. 2014; Finlator et al. 2011; Jaacks et al. 2012b; Lundgren 
et al. 2014). Most recently this has been examined by Salmon 
et al. (2014), who found that the star- formation rates of galax- 
ies from z = 3 to 6 are consistent with a power-law of the form 
\I/(t) = (t/r ) 7 (with 7 = 1.4 ± 0.1 and r = 92 ± 14 Myr). This 
analysis assumed that studying galaxies at a constant number 
density allows one to trace the progenitors and descendants 
of a galaxy population (e.g., van Dokkum et al. 2010; Leja 
et al. 2013), and their star- formation history was measured 
for a constant cumulative number density of 2 x 1 0 4 Mpc " 3 . 
Although the accuracy of this constant number density tech- 
nique was initially studied atz < 3, recent evidence shows that 
it likely works out to z ~ 8 (albeit it with a possible slight evo- 
lution of number density with redshift; Behroozi et al. 2013, 
Jaacks et al. in prep). 

Using our updated luminosity functions, we examine 
whether the star-formation histories at this earlier epoch are 
consistent with a similar functional form. Figure 20 shows 
the cumulative luminosity functions at z = 4 to 8 from our 
analysis. Using the Salmon et al. rising SFH, we can evolve 
our z = 7 cumulative luminosity function back in time to z = 8 
via: 

(9) 

where tp is the SFR, t : is the cosmic time elapsed since forma- 
tion to a given redshift, and using the Kennicutt (1998) con- 
version between Muv and SFR (with the updated coefficient 
of 1.25). The available data cannot constrain the formation 
redshift (zy), as it is degenerate with the star- formation his- 
tory exponent, thus we assume a value of zy = 1 8, which gives 
a close match between predicted and observed z = 8 cumu- 
lative luminosity functions. Figure 20 shows this predicted 
z = 8 cumulative luminosity function alongside our observed 
one. A very close match is seen at nearly all magnitudes. Our 
predicted z = 8 data points slightly under-predict the UV lu- 
minosity at M uv > -19. However, as discussed above, our 
constraints on the faint-end of the luminosity function at z= 8 
are tenuous at best. The agreement at the bright-end is excel- 
lent. While we did not correct for dust in this analysis, dust 
is highly unlikely to change these results (particularly at the 
bright end where we are interested), as bright/massive galax- 
ies at 4 < z < 7 all have similar UV slopes (Finkelstein et al. 
2012b; Bouwens et al. 2013). 

It is apparent when examining our cumulative luminosity 
functions in Figure 20 that this type of evolution will not work 
at all redshifts, as our luminosity functions are not uniformly 
spaced in magnitude (e.g., the z = 4 and 5, and z = 6 and 7 
cumulative luminosity functions are very close together). We 
examined one other redshift, evolving the observed z = 5 lu- 
minosity function to z = 6 . We find a decent match, though 
this under-predicts the bright-end, and over-predicts the faint- 
end. In any case, as we are most interested in extrapolating 
to z > 8, the fact that the predicted evolution works extremely 
well from z = 7 to 8 gives us confidence that extrapolating 
to higher redshifts is reasonable. This assumed evolution is 
stronger than that observed from z = 6 to z = 7. Had we as- 
sumed a SFH which matched the evolution from z = 6 to z = 7, 
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Fig. 20. — Left) The cumulative luminosity functions at z = 4 to 8 from our study. We evolve the observed luminosity functions to higher redshift, assuming 
that the Muv oc SFR, that the SFR rises with time as oc / I-4 (Salmon et al. 2014), and that galaxy progenitors and descendants share a common number density 
(e.g., Leja et al. 2013). We show two results: the predicted z = 6 luminosity function, evolved from z = 5 (blue squares), and the predicted z = 8 luminosity 
function evolved from z = 7 (yellow squares). Though there are small discrepancies at z = 6, the match between predicted and observed is excellent at z = 8. 
Right) Our differential luminosity functions at z = 7 and 8, with z = 9 data from the literature (triangles and circles). The red (gray) squares show our predicted 
z = 9 (10) luminosity function, continuing the evolution from z = 7 as shown in the left panel. The red line shows the z = 7 best-fit Schechter function, dimming 
My V to -19.95, providing an excellent match for our predicted z = 9 luminosity function. The dashed ride line shows this same function, but with a equal to 
the z = 8 value of -2.3. This predicted luminosity function shows excellent agreement with the observed faint-end at z = 9 when assuming the z = 8 value of a, 
but is significantly higher at the bright end compared to the published luminosity function of Oesch et al. (2013). However, the recent discovery of bright z = 10 
candidate galaxies by Oesch et al. (2014) (large gray squares) imply that bright galaxies are indeed present at this early epoch. However, we note that the z = 10 
number densities from Oesch et al. (2014) are actually more consistent with our z = 9 predictions than z = 10, thus clearly more work is needed to sort out the 
bright end at such early times. We conclude that when sufficient data exists for a large-volume survey for z = 9 galaxies, large numbers of bright galaxies will be 
discovered. 


we would have over-predicted the z = 8 LF. Our use of a SFH 
which matches the observed z= 7 to z = 8 evolution thus yields 
a conservatively low z = 9 predicted luminosity function. 

Given the relative paucity of observational information at 
z > 8, the fact that our assumed SFH matches the evolution 
from z = 7 to 8 makes it interesting to continue our study 
out to z = 9. Figure 20 shows the expected z = 9 luminosity 
function from our model, alongside our observed luminosity 
functions at z = 7 and 8. We calculated the expected z = 9 
luminosity function by again taking the z = 7 luminosity func- 
tion and evolving it out to z = 9 assuming the star-formation 
history discussed above (and Zf= 18 for all number densi- 
ties/magnitudes). As shown in this figure, our predicted z = 
9 luminosity function is consistent at the ~1 a level with all 
published data points from McLure et al. (2013) and Oesch 
et al. (2013). The insignificant under-prediction at the faint 
end is likely due to the fact that our analysis effectively keeps 
the faint-end slope fixed to the z = 7 value, while in reality 
it may become steeper. Figure 20 shows a Schechter func- 
tion which matches our predicted z = 9 luminosity function; 
we hold constant the value of a and ip* from z = 7, and find 
that M^v = -20.05 best matches our predicted data. This is 
more than a magnitude brighter than that reported by Oesch 
et al. (2013), yet still moderately consistent with the observed 
data points from both Oesch et al. (2013) and McLure et al. 
(2013) at the faint end. Using our predicted z = 9 luminos- 
ity function, we would expect to see 4.0 z = 9 galaxies at 
Muv < -20.3 ( H < 27) in a GOODS-sized field. Based on 
Poisson statistics alone, this is in mild tension with the zero 
galaxies at these magnitudes reported by Oesch et al. (2013). 
We also show results from this analysis when evolving z = 7 
to z = 10, which predict Mf} v = -19.7 (or ~1.5 Muv < -20.3 
galaxies per GOODS field). 


Recently, Oesch et al. (2014) have performed a new search 
for extremely distant galaxies, finding four bright z = 10 can- 
didates in the GOODS-N field and two new candidates from a 
re-analysis of the GOODS-S dataset. Although as mentioned 
above, these fields do not have deep F140W data, Oesch et al. 
(2014) used the < 3er detections of these galaxies in the ex- 
tremely shallow 800s 3D-HST (PI van Dokkum) F140W pre- 
imaging data to place these galaxies at z = 10. Even though 
these galaxies are only detected in one band with HST 6 , their 
presence is intriguing. Figure 20 shows the number densities 
of these sources from Oesch et al. (2014); there is excellent 
agreement with our predicted z = 9 evolution, though these 
data are much higher in abundance than our predicted z= 10 
luminosity function. Although these sources may be at z= 10 
rather than z = 9, if real, their presence confirms that bright 
galaxies are relatively abundant at z > 8.5. 

Finally, we examine the change in the integrated luminosity 
density at z = 9 with our proposed luminosity function com- 
pared to that from Oesch et al. (2012). A change only in My V 
makes little difference in the luminosity density, which is not 
surprising given the general shape of the luminosity function. 
However, as shown in Figure 20, our predicted z = 9 luminos- 
ity function does slightly under-predict the published number 
densities at the faint end. If we change the faint-end slope to 
-2.36 (to match our z = 8 value), we find a better agreement 
with the published faint z = 9 number densities. The lumi- 
nosity density with this steeper faint-end slope, including our 
brighter value of My V , is ~60% higher than that published 

6 Oesch et al. (2014) do claim detections of all four of the GOODS-N 
sources in the IRAC bands, but at least two of these sources are heavily 
blended with a nearby bright source, and it is unclear whether robust pho- 
tometry from a faint source can be recovered from such an environment. 


CANDELS: The Rest-Frame UV Luminosity Function at z = 4-8 
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in Oesch et al. (2012). Thus, it may be that the precipitous 
decline in the luminosity density (Oesch et al. 2013, 2014) 
(though see also, e.g., Ellis et al. 2013; Coe et al. 2013) may 
be less than previously thought (e.g., Ellis et al. 2013; Coe 
et al. 2013; Behroozi & Silk 2014). While these results are 
intriguing, we conclude that in order to robustly probe the 
bright end of the z = 9 luminosity function, we require a sig- 
nificantly increased searchable area with the correct filter set 
(allowing more than single-band detections) to discover these 
distant galaxies. Constraints on the full shape of the luminos- 
ity function in this distant epoch are crucial to design the most 
efficient surveys with JWST. 

10. CONCLUSIONS 

Combining the extremely deep data available in the HUDF 
with the still deep yet much wider data available from CAN- 
DELS in the GOODS-South and GOODS-North fields allows 
robust samples of galaxies to be discovered across a large dy- 
namic range of UV luminosity at z = 4, 5, 6, 7 and 8. Using a 
robust photometric redshift selection technique, we have dis- 
covered a sample of nearly 7500 galaxies at3.5<z<8.5 
over five orders of UV magnitude, and over a volume of 0.6- 
1.2 x 10 6 Mpc 3 . We discovered a large number of bright 
(A/uv < -21) galaxies at these redshifts, in excess of predic- 
tions based on previous estimates of the luminosity functions 
atz > 6. 

• Our sample selection performs very well when com- 
paring to available spectroscopic redshifts. We perform 
various tests to estimate the contamination rate, which 
we find at worst to be <15%, and more likely to be 
<5-10%. This is consistent with contamination esti- 
mates based on the colors of the most likely contami- 
nants, dusty star-forming galaxies at z ~ 2. Although 
the GOODS fields are only two of five CANDELS 
fields, the remaining three fields contain relatively shal- 
low 7-band data, which can result in increased sample 
contamination, as well as a reduced ability to separate 
galaxies into z = 6, 7 and 8 samples. 

• Our large volume probed allows to us make a robust 
determination of the amplitude and shape of the bright- 
end of the luminosity function, which can be used as 
a crucial probe of the physics dominating galaxy evo- 
lution. We used a MCMC technique to estimate the 
luminosity function, to better characterize the uncer- 
tainties, both on the step-wise luminosity function, as 
well as on the parameters of the Schechter functional 
form. Our results agree with previous studies at the 
faint end, but deviate from some previous studies at the 
bright end, where our data allow us to better constrain 
the abundance of rare, bright galaxies. We find results 
consistent with a non-evolving characteristic magnitude 
(My V « -21). Both the faint-end slope (a) and the nor- 
malization (<£>*) do significantly evolve with increasing 
redshift, to steeper and lower values, respectively. This 
is in contrast to previous results, which determined that 
the evolution of the luminosity function was primarily 
in luminosity. 

• We explored whether a Schechter functional form is 
required by the data, or whether a single (or double) 
power-law is a better fit for our luminosity functions; a 
single power-law form of the luminosity function may 


be expected at very high-redshift, when dust may not 
be present, and/or feedback due to AGN activity is no 
longer sufficient to suppress star-formation in the most 
massive galaxies. At z = 6 and 7, a Schechter (or double 
power-law) is required to fit the bright end. Flowever, 
at z = 8, a single power-law provides an equally good 
fit to the data. Although larger volumes will need to 
be probed to improve the estimates of the abundances 
of bright z = 8 galaxies, if a power law is preferred, 
it could imply that we may be observing the era when 
feedback stops affecting massive galaxies. Comparing 
to semi-analytical models, we find that the evolution in 
our luminosity function can be explained by a chang- 
ing impact of dust attenuation with redshift. In a fu- 
ture work we will explore whether this is a unique con- 
straint, or whether a combination of feedback and dust 
changes can reproduce the observations. 

• We measure the evolution of the cosmic star-formation 
rate density by integrating our observed luminosity 
functions to the observational limit of Myy = -17, and 
correcting for dust attenuation. We find that the cos- 
mic SFR density evolves as (l+z)~ 4 - 3±a5 atz > 4. This 
smoothly declining function with increasing redshift is 
consistent with published estimates of the SFR density 
at z > 9. 

• We investigate the constraints on the contribution of 
galaxies to reionization by integrating our luminosity 
functions down to M\jy = -13. We find that our fiducial 
results (assuming C/ f esc = 23, which does not violate 
Lyo forest constraints at z < 6) are consistent with a 
reionization history that begins at z > 10, and completes 
atz « 6, with a midpoint at 6.7 < z < 9.4. Flowever, the 
uncertainties, particularly at z > 7 are high, due to the 
relatively high uncertainty in the faint-end slope, such 
that our observations are consistent with an IGM at z = 8 
that is anywhere from completely ionized, to 85% neu- 
tral. 

• The presence of bright galaxies at z = 6 - 8 has interest- 
ing implications for the luminosity functions at higher 
redshift. We used empirically derived star-formation 
histories to evolve ourz = 7 luminosity function back to 
z = 9, and predict that ~4 bright (Mjjy < -20.3) galax- 
ies should be detectable per GOODS-sized field. This 
is contrary to initial observational results, which, us- 
ing single-band detections found no bright z = 9 galax- 
ies, though consistent with emerging results that some 
bright galaxies may exist at z = 10. Future wider-area 
studies with two-band detections will provide a more 
robust estimate of the bright-end of the z= 9 luminosity 
function. 

This study highlights the power of combining deep and 
wide-area studies to probe galaxy populations at very high 
redshifts, a topic that will remain highly active through the 
advent of JWST. These results leave us with a variety of ques- 
tions. What is responsible for the apparent abundance of 
bright galaxies at z > 6? Is this tied in with a reduction of 
feedback, or is some other physical process in play? Does 
this trend continue out to higher redshifts, or does the lumi- 
nosity density fall off dramatically at z > 8, as has been pro- 
posed? Although these issues are inherently intertwined, we 
can make progress on these issues now with future wide-area 
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HST surveys. This will allow us the most complete view of 
the high-redshift universe by the end of this decade, allowing 
us to make full use of JWST. 
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